Endothelium and Subendothelial Matrix Mechanics Modulate Cancer Cell Transendothelial Migration

Abstract Cancer cell extravasation, a key step in the metastatic cascade, involves cancer cell arrest on the endothelium, transendothelial migration (TEM), followed by the invasion into the subendothelial extracellular matrix (ECM) of distant tissues. While cancer research has mostly focused on the biomechanical interactions between tumor cells (TCs) and ECM, particularly at the primary tumor site, very little is known about the mechanical properties of endothelial cells and the subendothelial ECM and how they contribute to the extravasation process. Here, an integrated experimental and theoretical framework is developed to investigate the mechanical crosstalk between TCs, endothelium and subendothelial ECM during in vitro cancer cell extravasation. It is found that cancer cell actin‐rich protrusions generate complex push–pull forces to initiate and drive TEM, while transmigration success also relies on the forces generated by the endothelium. Consequently, mechanical properties of the subendothelial ECM and endothelial actomyosin contractility that mediate the endothelial forces also impact the endothelium's resistance to cancer cell transmigration. These results indicate that mechanical features of distant tissues, including force interactions between the endothelium and the subendothelial ECM, are key determinants of metastatic organotropism.


Supplementary note 1-Pore size and fiber thickness of collagen matrix
To accurately quantify the fiber thickness and pore size we imaged the dehydrated samples (using graded ethanol) and employed a scanning electron microscope (Zeiss-Gemini, Germany). Collagen is the main contributor to the stiffness and strength of soft connective tissues, such as the cornea, blood vessels, ligament and tendons [1] , in which the uniaxial alignment, high collagen density, hierarchical fibril architecture, and degree of crosslinking enables these tissues to resist high tensile loads [2] . However, collagen hydrogels represent relatively poor mechanical properties in the absence of covalent cross-linking [3] and increasing collagen content does not seem to improve these properties tremendously. Therefore, applying some modifications, such as using crosslinkers, seems to be more physiologically relevant than elevating collagen concentration. Furthermore, in contrast to some other hydrogels present in native ECM, collagen hydrogels do not have any growth factor (GF)-specific binding sites and often cannot sequester physiological amounts of the protein [3] and again employing chemical cross-linkers would provide a physiologically relevant strategy to improve the hydrogel GFbinding properties.
Such limitations have triggered many studies on collagen gels to provide a more physiologically relevant platform disassociating the effects of hydrogel stiffness from other factors. For instance, to de-couple matrix stiffness from matrix density and structure in collagen gels, Mason et al. [4] used nonenzymatic glycation of the collagen in solution, prior to polymerization to increase the compressive modulus of the hydrogel without significant changes to its architecture. They demonstrated that matrix stiffness alone can promote spreading of endothelial cells within the 3D matrix and consequently, number and length of the angiogenic sprouts enhanced. Additionally, carbodiimides [5] , polyethylene glycol [6] , glutaraldehyde [7] , isocyanate [8] , glycation with glucose [9] or ribose [4] , and genipin [10] have been used as chemical cross-linkers to improve mechanical properties of collagen hydrogel. An interesting avenue to explore in future studies would be the use of crosslinkers to tune mechanical properties of the matrix without altering its structural architecture to provide a better picture of the mechanical behaviour of EC / TC during transmigration.
Gelation temperature has a significant effect on polymerization kinetics that will largely determine the mechanical properties of collagen hydrogels. Fibrillogenesis occurs faster at higher temperatures due to accelerated nucleation and lateral aggregation of collagen molecules [3] . Shannon et al. [11] showed that for 3 mg/ml collagen I concentration, increasing gelation temperature from 18 °C to 37 °C leads to a reduction in polymerization time. Similar observations were reported by other researchers [12,13] . Furthermore, our SEM images show that reducing the temperature results in the formation of thicker fibers. Other researchers also observed a similar pattern using SEM [12,14,15] , confocal microscopy [16] , confocal reflectance microscopy [17] , multiphoton microscopy [15] , and second harmonic generation microscopy [18] . It has also been shown that fibers polymerized at higher temperatures are shorter, randomly aligned, and less bundled [19] . Lower temperatures are thought to limit nucleation of new fibers via decreasing entropy, which promotes thickening and elongation of already existing fibers, and forms networks, which are often more heterogeneous [3,20] . Additionally, in agreement with previous studies [15][16][17]19,21] , our SEM images show that pore size increases with decreasing gelation temperatures for the same level of collagen content (Fig 1f, S1, S2). Additionally, a higher permeability for the porous gel would provide additional evidence indicating its larger pore size (Fig S3).
However, regarding mechanical properties, there is no agreement among researchers; some studies have shown that higher gelation temperature elevates the stiffness [21] , while other researchers have shown that greater compressive or shear moduli are produced in gels at lower temperatures. Achilli and Mantovani [22] reported that for pH=7, increasing the temperature increased stiffness, whereas a reverse pattern was observed for pH=10. These discrepancies are thought to depend on pH and collagen concentration, and possibly on the mechanical testing setup [3] . Fig S1) Scanning electron microscopy image of the control, high collagen content, and porous gels (Scale = 1 um) Fig S2) Quantification of pore size and pore area for hydrogels. a) SEM image from the sample. b) binarisation of the image using the "Trainable Weka Segmentation" tool in ImageJ. Only fibers located on top surface (red) were defined as region 1 while corresponding pores (green) were selected as region 0. c) Binarized image. d) Pores were labelled, and their areas were calculated using the "Analyze Particle" tool in ImageJ. e) The average area of pores (mean ± s.e.m., n=3 samples, in each sample, pore size evaluated at 3 randomly selected regions; data points represent the average pore area calculated for each dish; *p < 0.05, **p<0.01, ***p<0.001, Student's t-test). f) The biggest circle that could be inscribed in each pore (green circles) was found using MATLAB. Scale =1 μm.

Device Fabrication for permeability measurement
Acrylic moulds were fabricated by cutting (Epilog laser cutter) our designs through 1 mm thickness acrylic sheets. Next, the laser-cut acrylic geometries were bonded to larger acrylic bases using acrylic cement (Weld On 16 Clear, Multibond Solutions). These moulds were then taped to the bottom of Petri dishes to fabricate the final Polydimethylsiloxane (PDMS) devices. PDMS and crosslinker (10:1; Sylgard-184; Dow Corning) were mixed, degassed, and poured over the mould up to a height of 4 mm. After curing for 2 hours at 90° C, the PDMS was cut and peeled from the mould. Ports at the ends of the channel were punched out using 2 mm and 4 mm biopsy punches (Miltex). The device was cleaned with tape to remove debris. Then the bottom of the device and a glass coverslip were treated with plasma for 60 s and pressed firmly together to form an irreversible bond. To restore hydrophobicity of surfaces and create irreversible bonding, the final devices were placed inside an 80°C oven overnight. Next, the collagen solution was injected to the collagen channel (length = 25 mm, width = 3 mm, height = 1 mm) from the smaller port, followed by hydrogel incubation at 37°C or 20°C. Finally, to measure the permeability, 50 μL dye was added to the larger port and the devices were imaged every 15 min. Fig S3) Quantification of the permeability of a dye in control and porous gels. a) Time lapse images showing the penetration of the blue dye along a gel channel made in PDMS. Scale = 1 cm. b) Distance travelled by the dye (x) vs time. c) Comparison on penetration length after 2 hr for control and porous gels (Data points represent the data obtained from each experiment, n=3 independent experiments, mean ± s.e.m.; *p < 0.05, **p<0.01, ***p<0.001, Student's t-test).

Fig S4)
Elastic modulus the collagen gels measured by atomic force microscopy indentations on the top surface of the gels. Collagen gels were prepared in a glass bottom dish and were incubated in 37 °C (control and HCC) or 20 °C (Porous) for two days, then, elastic modulus was measured using AFM without introduction of any cells to the dish (mean ± s.e.m., n =3 dishes. In each dish at least 40 measurements were carried out. Data points represent average elastic modulus of each gel; *p < 0.05, **p<0.01, and ***p<0.001, Student's t-test).

Supplementary note 2-Quantification of stiffness of the subendothelial matrix
In this note, we explain how the Elastic modulus of the subendothelial collagen matrix was calculated as a function of depth. The stiffness of the collagen gel was assumed to be uniform before seeding of the endothelial cells (ECs) on top. However, once formed, the endothelial monolayer remodeled the underlying collagen gel by applying shear and compressive forces. This led to a denser and stiffer layer of collagen fibers at the top of the subendothelial matrix. The next step was to characterize the stiffness of collagen at increasing depths through the gel. To achieve this, we used the following information: (1) Einitial: the stiffness of the collagen gel was measured before adding the ECs, using atomic force microscopy (AFM, Fig S5).
(2) Eapparent: the apparent stiffness of the collagen was quantified using AFM after decellularization of the EC monolayer using Triton X. After washing away the monolayer of ECs, the collagen gel swelled and its thickness increased, but it remained thinner than its initial thickness before any EC remodeling had occurred. Therefore, the gel was still non-uniform in depth, and we expected it to be stiffer near its top surface. Thus, we used the word "apparent" to emphasize that the gel does not have a single stiffness value under this condition. Eapparent represents an average stiffness of the layers near the top surface and since E decreases with increasing depth, Eapparent should be lower than the Elastic modulus at the top surface of the gel (Etop).
(3) confocal images from collagen fibers before and after washing off the EC monolayer. This enabled us to quantify "mesh size" (z) [23] of the gel at each z-stack plane of focus. The mesh size is an index that represents the average distance between collagen fibers. From (3), we can calculate the mesh size at each z-plane of focus. We then needed to devise a method to correlate the Young's modulus with the mesh size to be able to estimate the Young's modulus at each plane of focus. As suggested in previous studies [23] , the Young's modulus is inversely proportional to a power of mesh size: in which A and n are two constants. To find these two parameters, we needed to know the values of the Young's modulus and the mesh size at two focus planes, such as the top and bottom of the gel. Since the bottom is far away from the EC monolayer, it would not be affected by the EC monolayer and the Young's modulus would remain equal to the initial value (Ebottom=Einitial). However, determining Etop is not quite as straightforward and requires trial and error. In the following section, we provide a step-by-step procedure to calculate the Young's modulus for the porous gel.

Determining mesh size at increasing depths
As mentioned above, mesh size is an index representing the void size of the hydrogel. To find the mesh size (z), we used the protocol proposed by Yang et al. [23] . Briefly, z-stack confocal images at each plane of focus were binarized after being segmented using a thresholding algorithm in ImageJ. Then, the horizontal and vertical distances between collagen fibers were measured using a MATLAB code. The mesh size can be obtained from the exponential fit to the void size histogram (Fig S6). The horizontal and vertical void sizes show similar distributions. This provides further evidence that the collagen gel can be considered as an isotropic material at each focal plane. (− . ( ))). c) By repeating this procedure, the mesh size can be calculated at increasing depths through the gel (Data are mean ± s.e.m., n= 8 samples).
Finally, using an exponential fit, the mesh size can be written as a function of depth. For the porous gel:

Finding the Young's modulus at the top and the exponent n
For the porous gel the initial and apparent Young's moduli were measured using AFM (Einitial=Ebottom=69 Pa, Eapparent = 100 Pa, Fig S7a). Additionally, the mesh sizes were calculated in the previous step at all focal planes, in particular, at the top and bottom of the gel (ztop=10, zbottom=16.7). Our next aim was to estimate the Young's modulus at the top of the gel (Etop) by trial and error. To achieve this, we carried out the following steps: 2) Calculate the exponent n (from Eq. S1) for each guesstimate value. Since the Young's moduli and mesh sizes are known at the top and bottom of the gel, n can be found using the following equation for each case: Using Eq. S3, the exponent n was found to be 1.5, 1.0, and 0.6 for Etop=190, 135, and 100 Pa, respectively. 3) Combining Eq. S1, S2, and the value found for n, the Young's modulus can be determined as a function of depth (z) for each guess (Fig S7c): We employed an axisymmetric model in which the young's modulus of the gel was defined as a function of depth (given by Eq. S4), while the indentor was rigid (Fig S7d).
5) The force indentation curves obtained from AFM experiments were compared to that from the finite simulations for each guess (Fig S7e). The Young's modulus at the top of the gel and the corresponding exponent n would be adopted from the best fit (For the porous gel, Etop=135 Pa and n=1.0).

Quantifying the depth dependent Young's modulus before washing the EC monolayer
Having quantified the exponent n, from the decellularized gel, we next calculated the Young's modulus of the collagen matrix before decellularization and washing off the EC monolayer. To this end, the following steps were carried out: 1) Binarizing the confocal images taken from collagen fibers (before washing off the EC monolayer) and quantifying the mesh size at each focal plane (Fig S8a, b).
2) Calculation of the Young's modulus as a function of depth using the following equation (Fig S8c): Such procedures were repeated for control and HCC gels to characterize the mechanical properties of these gels as a function of depth.

Immunofluorescence staining for MLC 2 and p-MLC2
Cells were fixed using 4% PFA at room temperature for 15 min. Cells were then washed in PBS, permeabilized for 10 min with 0.3% Triton-X in PBS and blocked for 1 h in blocking buffer (10% normal goat serum and 0.1% Triton-X in PBS). The primary antibodies were then applied overnight at 4 °C with gentle shaking. The following primary antibodies were used: Myosin Light Chain 2 (MLC2, Insight Biotechnology, 10906-1-AP) at a concentration 1:50 and Phospho-Myosin Light Chain 2 (pMLC2, Cell Signaling Technology, 3671S) at a concentration 1:50. Cells were then washed and incubated with secondary antibodies at room temperature for 2 h. The following secondary antibody was used: goat anti-rabbit (Cell Signaling Technology, 4412S) at a concentration of 1:500. Finally, cells were washed and stained with Hoechst 33342 for 10 min at room temperature. Samples were imaged using an upright confocal microscope (Zeiss Airyscan LSM 980). To find the area below the cell, we used the CD31 (yellow) channel to observe the boundaries of each cell. Then, after segmentation, the boundaries were skeletonized, and tractions were integrated inside each boundary. Cells in which their area was too big or too small were dismissed. C) Validation of the model. To quantify errors involving our force quantification technique, we quantified the forces generated by an isolated cell. Since the cell was stationary, the net force being applied to it should be zero. As there is no other cell in the vicinity of the cell, the tugging force is zero and the integral of traction forces should also be zero. Non-zero values show the level of error in our method which is 0.28 ±0.09 nN. The error is less than 16% of the minimum calculated forces (for EC-RhoA case) and less than 6.5% of forces measured for the control case (p-Value =1.3 e-6 and 4.7e-11, respectively; mean ± s.e.m, n ≥ 8 samples, ***p<0.001, Student's t-test)

Supplementary note 4-Details of the computer simulations
To simulate the mechanical interaction between the EC monolayer and ECM, we developed a computational model using the commercial finite element software ABAQUS in which a layer on ten cells was placed on top of a substrate. For the sake of simplicity, each cell was modelled as a rectangle while an ellipse was used to describe its nucleus. Typical values were chosen as the dimensions of cells (20 µm x 5 µm) and their nucleus (16 µm x 3 µm as the major and minor axes, respectively). Also, to avoid boundary effects, dimensions of the gel were selected to be sufficiently large (200 µm x 100 µm). Mesh size was also adopted after mesh sensitivity analysis. Additionally, the cells and the underlying gel were constrained in the lateral direction on the left, while the bottom of the gel was fixed in the vertical direction (Fig S13a). To simplify the model and reduce the computational time, a plane strain formulation was used to solve the boundary value problem. A linear elastic material model was employed to describe the mechanical properties of the gel and EC nuclei. A typical Young's modulus (500 Pa) and Poisson's ratio (0.4) were selected as the material properties of the EC nucleus, as these parameters were shown to have a negligible effect of the computational results. For the collagen substrate, Poisson's ratio was 0.2 [24] while a range of values (100 -500 Pa), that covers the stiffness of the collagen gels used in this study, was selected as its Young's modulus. To describe the mechanical properties of the cytoskeleton, a model developed by Shenoy et al. [25] was used, in which cells are considered as materials with passive and active properties. As a passive material, cells would deform under external forces while the active element enables them to apply contractile forces to their surrounding microenvironment. To consider the passive behaviour, the model uses two parameters (bulk modulus, K, and shear modulus, ) while to accommodate the active element, the model requires three more parameters (motor density in the quiescent state, 0 , chemical stiffness, , and chemomechanical feedback parameter, ). 0 shows the contractile ability of an isolated cell in the absence of ECM or any other cell, controls the effect of molecular mechanisms that regulate the engagement of motors, while is responsible for the feedbacks affecting the stress-dependent signalling pathways. As porosity of the ECM regulates actin polymerization and RhoA signalling, its effect was taken into account through changing the feedback parameter . RhoA activity and cell generated force are less, when we use a porous substrate, thus a lower a should be used for more porous matrix. All the cytoskeletal material parameters were taken from Shenoy et al. [25] , except , as a slight increase in its value led to us obtaining more accurate results. Material parameters are also shown in (Fig S13a). Fig S13. Computer simulation details. a) Geometry, boundary condition, material models, and material parameters used for computer simulation in ABAQUS. b) Quantifying horizontal and vertical stiffness (i.e., force required to generated unit displacement) Formulation of the model for 3D condition has been provided previously by Shenoy et al. [25] , briefly: in which, is the stress tensor, is the Kronecker delta and represents the contractility tensor. = 11 + 22 + 33 and = 11 + 22 + 33 where is the strain tensor.
For the plane strain case, the formulation can be simplified as below: where 12 = 2 12 . The constitutive model (i.e., SI Eq. 2 and 3) were implemented into ABAQUS using a UMAT developed in Fortran.
To obtain the horizontal and vertical stiffnesses, i.e., the forces required to generate unit displacements either in lateral or vertical directions, we modified the geometry slightly, while the other dimensions, material models and properties remained unchanged. We created a 5 µm gap on the left-hand side of the monolayer; this allowed us to apply a 1 µm horizontal displacement to the left of the monolayer and measure the concomitant reaction forces. Furthermore, we applied a 1 µm vertical displacement to the gap and quantified the vertical stiffness of the matrix (Fig S13b).

Supplementary Note 5-Mechanics of cancer extravasation
The aim of this note is to illustrate the traction distribution on a representative cell for each condition. Also, average tractions and the duration were quantified for each stage. Impact of stiffness of the collagen substrate (E) and its porosity (through feedback parameter a) on the lateral stiffness of EC monolayer (i.e., the force that needs be applied to ECs to generate unit lateral displacement in the monolayer) and vertical stiffness (i.e., the force that needs be applied to generate unit vertical displacement in ECM) obtained from computer simulations. Points A, B, and C denote conditions corresponding to control, high collagen content (HCC), and porous cases, respectively. The same feedback parameter a (4.64E-3 Pa -1 ) was used to simulate the behaviour of collagen gel in control (point A) and HCC (point B) conditions as the mesh size is not significantly different for these cases. However, since for the porous condition, the mesh size is larger than that for the control, a smaller a (4.2E-3 Pa -1 ) was employed for that case. Additionally, the Young's modulus of collagen gel for the control and porous gel was in same range (164 Pa), while a higher value (460 Pa) was used for the HCC condition. The contour plots show that the lateral stiffness of the monolayer for the HCC case (~1700 µN/m) is larger than that for the control case (~1100 µN.m) while the monolayer formed on the porous gel has the minimum lateral stiffness (~650 µN/m). Similarly, the vertical stiffness of the gel in the HCC case (~525 µN/m) is more than that for the control case (~210 µN/m) and the minimum vertical stiffness belong to the porous case (~190 µN/m).

Supplementary Note 6-Impact of EC density
For all experiments in this study, EC monolayers were formed by introducing 350,000 HUVECs to the MatTek dish. We chose this EC density after optimizing the culturing protocol for ECs at different densities (2.e5, 3.e5, 3.5e5, and 4.e5) and investigating the confluency and integrity of the monolayer after 48 hrs. Using standard T75 flasks, we estimated that the initial cell density should be ~3. -4.e5 per dish, since ~3e6 cells provide a confluent monolayer in the flask, and the area of our dish is ~10.2 cm 2 (~7 times less than a T75 flask). Considering the proliferation rate of ECs, using an initial density of ~3. -4.e5 cells should result in a seamless confluent monolayer after 48 hrs. The results showed that with 3.e5 cells, there were some non-confluent regions in the dish, while for the density of 5.e5, we found some overconfluent regions (Fig S22)." Fig S22) Effect of endothelial cell density on confluency of the monolayer. The monolayer with 2e5 and 3e5 initial cell densities, was not fully confluent, while for 5e5 cells, some regions were overconfluent. Using a cell density of 3.5e5, a seamless confluent EC monolayer was formed after 48 hr of culture. Scale = 200 μm.
To investigate the impact of cell density on transmigration potential, we formed EC monolayers with two different initial densities: 3.5e5 and 7e5. The results revealed that for the higher EC density, transmigration rate and efficiency decreased significantly, implying that it is more difficult for TCs to transmigrate through a compacter monolayer (Fig S23). Fig S23) Effect of EC density on a) transmigration rate and b) efficiency (mean ± s.e.m., n=3 samples (dishes); data points represent the average transmigration rate / efficiency for each dish. *p < 0.05, Student's t-test).

Supplementary note 7 -Verification of the force microscopy method
The aim of this note is to evaluate the accuracy of the method used to quantify the cell-generated tractions. The protocol included the following steps: 1-Calculating the Young's modulus of the subendothelial matrix as a function of depth, (supplementary note 2).
2-Taking time-lapse confocal images of the collagen fibers of the subendothelial matrix. At least two time-points, t and t+Dt, are required.
3-Measuring the displacement field in the matrix between t and t+Dt using the code developed by Bar-Kochaba et al. [26] .
4-Calculating the increment of cell-generated tractions between t and t+Dt, by implementing the data obtained from the two previous steps into a commercial finite element software (ABAQUS).
Here, we used a collagen gel (6 mg/ml, incubated at 37 o C for 30 min) and quantified its Young's modulus as a function of depth after washing the EC monolayer, using the method explained in Supplementary note 2. The Young's modulus and mesh size were 300 Pa and 1.67 at the bottom of the gel, respectively. Considering the mesh size at the top (x=1.34), the Young's modulus was calculated at the top (Etop = 700 Pa, Fig S24). To verify the accuracy of the method, we must apply a known force to the matrix, then follow the above-mentioned steps to estimate the force and compare the result with the applied force. To this end, we employed AFM experiments to indent the gel shown in Fig S25, using a spherical tip glued to the end of cantilever. The sample stage of the AFM was mounted on an inverted confocal microscope, allowing us to image fluorescently labelled collagen fibers at the same time. We applied a 35 nN force to the gel which caused 8.2 µm indentation depth. Then, the following steps were carried out to reconstruct the force: 1) Confocal images of the collagen fibers were captured (Fig S25a), 2) Calculating the displacements using the FIDVC code (Fig S25b), 3) Performing the finite element simulation and estimating the stress distribution on the AFM tip (Fig S25c) and, 4) Integrating the stresses to measure the total force applied by the AFM tip to the matrix and compare the reconstructed force with AFM measurements (Fig S25d). Comparison of the reconstructed normal (28.1 nN) and shear (5.9 nN) forces with the applied forces (35 nN and 0 nN, respectively) shows that the error in prediction of forces could be up to 20%, which is within an acceptable range.