Cell Shape and Durotaxis Explained from Cell-Extracellular Matrix Forces and Focal Adhesion Dynamics

Summary Many cells are small and rounded on soft extracellular matrices (ECM), elongated on stiffer ECMs, and flattened on hard ECMs. Cells also migrate up stiffness gradients (durotaxis). Using a hybrid cellular Potts and finite-element model extended with ODE-based models of focal adhesion (FA) turnover, we show that the full range of cell shape and durotaxis can be explained in unison from dynamics of FAs, in contrast to previous mathematical models. In our 2D cell-shape model, FAs grow due to cell traction forces. Forces develop faster on stiff ECMs, causing FAs to stabilize and, consequently, cells to spread on stiff ECMs. If ECM stress further stabilizes FAs, cells elongate on substrates of intermediate stiffness. We show that durotaxis follows from the same set of assumptions. Our model contributes to the understanding of the basic responses of cells to ECM stiffness, paving the way for future modeling of more complex cell-ECM interactions.


INTRODUCTION
Mechanical interactions between cells and the extracellular matrix (ECM) are crucial for the formation and function of tissues. By sensing and responding to physical forces in the ECM, cells change their shape and migrate to other locations. The shape of a wide range of mammalian cell types depends on the stiffness of the ECM. In vitro, cells cultured on soft, two-dimensional ECM substrates become relatively small and round, whereas on stiffer ECMs the cells assume elongated shapes. On highly rigid substrates such as glass, cells spread out and flatten. This behavior has been observed for a wide range of cell types, including endothelial cells (Califano and Reinhart-King, 2010), fibroblasts (Ghibaudo et al., 2008;Pelham and Wang, 1997;Prager-Khoutorsky et al., 2011), smooth muscle cells , and osteogenic cells (Mullen et al., 2015). Cells tend to migrate toward stiffer parts of the ECM, a phenomenon known as durotaxis. Such behavior also occurs for a wide range of mammalian cell types, including fibroblasts (Lo et al., 2000), vascular smooth muscle cells (Isenberg et al., 2009), mesenchymal stem cells (Vincent et al., 2013), and neurons and glioma cells (Bangasser et al., 2017). Here, we show that a single model suffices to explain (a) increased cell area on stiffer substrates, (b) cell elongation on substrates of intermediate stiffness, and (c) durotaxis.
It is still not completely understood what molecular mechanisms regulate such cellular response to ECM stiffness (Jansen et al., 2015). Cells are able to sense matrix stiffness through focal adhesions (FAs). FAs are multimolecular complexes consisting of integrin molecules that mediate cell-ECM binding and force transmission and an estimated further 100 to 200 protein species that strongly or more loosely associate with FAs (Bershadsky et al., 2003a;Winograd-Katz et al., 2014;Zaidel-Bar et al., 2007). Among these are vinculin and talin, which bind integrin to actin stress fibers.
FAs dynamically assemble and disassemble, at a rate that is regulated biochemically and mechanically. The disassembly rate is highest on soft ECMs and is lower on stiffer ECMs (Pelham and Wang, 1997) such that FAs stabilize more easily on stiffer ECMs. This mechanosensitivity of FA dynamics is regulated by talin and p130Cas, proteins that change conformation in response to mechanical force (Bershadsky et al., 2003a;Jansen et al., 2015). For instance, stretching the structural protein talin reveals vinculin-binding sites, allowing additional vinculin to bind to FAs (Rio et al., 2009) and stabilize the FA .
Interestingly, manipulations of FA formation affect both cell spreading and cell motility. Systematic knockdown of protein-tyrosine kinases involved in FA formation and traction force development reduces the ability of fibroblasts to elongate in a substrate-stiffness-dependent manner (Prager-Khoutorsky et al., 2011). FA assembly and disassembly has also been associated with cell migration and cell orientation (Broussard et al., 2008;Chen et al., 2013;Kim et al., 2013;Plotnikov et al., 2012) in response to the ECM. Hence, the mechanosensitive growth of FAs is key to our understanding of cellular responses to ECM stiffness, but how FAs and cytoskeletal force generation work together to regulate cell spreading, cell shape, and durotaxis is still to be elucidated.
Previous mathematical models for mechanosensitivity have proposed explanations for cell spreading, cell elongation, or durotaxis, but did not explain the three phenomena in unison. A central idea, dynamic mechanoreciprocity (Bissell et al., 1982;van Helvert et al., 2018), posits that properties of the ECM affect the mechanical and/or biochemical properties of the cell, that, in turn, change the local environment. For example, Ni et al. (Ni and Chiang, 2007) assumed that cells pull more strongly due to mechanical resistance of the matrix (Ni and Chiang, 2007), resulting in spreading of the cell; Shenoy et al. (Shenoy et al., 2015) proposed that stresses in the ECM can polarize the cytoskeleton. Cell spreading is also the main process addressed in several kinetic models of FAs (Deshpande et al., 2008;Ronan et al., 2014;Stolarska and Rammohan, 2017;Vernerey and Farsad, 2014). In these models it is assumed that the greater stress that FAs experience on stiffer substrates accelerate the formation of stress fibers and enable yet greater cell forces to be applied. The resulting positive feedback loop drives cells to spread out on stiffer substrates. More recently, such frameworks have been expanded with thermodynamic considerations (McEvoy et al., 2018;Shishvan et al., 2018). The prevalent idea of durotaxis is that increased stabilization of FAs at the stiff side of the cell leads to increased traction force (Dokukina and Gracheva, 2010), modified stress fiber dynamics (Harland et al., 2011;Lazopoulos KA, 2008), cell speed, and/or polarization (Kim et al., 2018;Malik and Gerlee, 2019;Novikova et al., 2017;Stefanoni et al., 2011;Yu et al., 2017), or other cell dynamics (Allena et al., 2016;Aubry et al., 2015;Shenoy et al., 2015).
Another explanation for cell anisotropy and durotaxis is based on the motor-clutch hypothesis (Chan and Odde, 2008): on stiff substrates, the cellular traction forces quickly break new adhesions, leading to a slip regime. On soft substrates, the cells can pull back while building up traction forces, gradually increasing load on the adhesions until they fail (Chan and Odde, 2008). Coupling multiple motor-clutch units, cell anisotropy emerges with an optimum that depends on the traction force and traction rate as well as on the number of motor-clutch units (Bangasser et al., 2017). Also, durotaxis emerges at a single cell (Bangasser et al., 2017) or collective (Sunyer et al., 2016) cell level.
Altogether, previous models suggest that mechanosensing is regulated by positive feedback between substrate stiffness, traction force, and/or FA stability. These models included many feedback mechanisms, with detailed dynamics of cytoskeleton/F-actin (Deshpande et al., 2008;Vernerey and Farsad, 2014), adhesion, and/or thermodynamics principles (McEvoy et al., 2018;Shishvan et al., 2018), resulting in a multiplicity of variables and parameters. On the plus side, such detailed models can be used to fit specific experimental data well and to incorporate facts from the experimental literature. On the minus side it remains unclear how the different mechanisms work together, which ones are essential or redundant, and what their respective roles are in distinct cell behavior such as elongation versus uniform spreading. Can durotaxis emerge from the same set of assumptions, or are additional mechanisms required?
In many models, it is assumed that cells become more contractile on stiffer substrates. However, it has been shown in fibroblasts that contractile forces do not depend on substrate stiffness (Feld et al., 2020;Freyman et al., 2002). Here we show that mechanosensitive focal-adhesion turnover, together with cell contractility and substrate adhesion suffice to explain all three ECM-stiffness-dependent cell behaviors. Our model suggests that cells distinguish between soft and stiff substrates not by adjusting the force magnitude, but rather the rate of force buildup: forces on FAs increase much faster on stiff substrates than on soft ones.
The model is based on the following assumptions: (1) The ECM is a nonlinear elastic material, such as a polyacrylamide gel.
(2) FAs are discrete clusters of integrin-ECM bonds. iScience Article the FAs from a pool of available free integrin bonds on the cell membrane. We assume a constant rate once the first adhesion bonds bring the cell membrane close enough to the ECM Sun et al., 2009. (4) The unbinding rate is suppressed by the tension in the FA due to pulling of stress fibers . (5) On soft ECMs it takes more time for the tension in the FA to build up to its maximal value, than on stiff ECMs . (6) Larger FAs detach less easily from the ECM. (7) Planar stress reinforces FAs due to recruitment of stabilizing proteins such as vinculin.
As we will show, the model predicts that the range of stiffness on which cells elongate depends on the velocity of myosin motor proteins. We find that simulated cells exhibit durotaxis. Consistent with experimental observations, the durotaxis speed increases with the slope of the stiffness gradient (Isenberg et al., 2009;Vincent et al., 2013). Because the key to our proposed model is mechanosensitivity of the FAs, it is likely generalizable to nonlinearly elastic, inhomogeneous, and fibrous natural matrices.

Modeling Approach and Details
We have extended our hybrid Cellular Potts-Finite Element framework  to include dynamic descriptions of FAs. Figure 1 gives a flowchart of the model, showing the feedback between the cell, the FAs, and the elastic substrate it adheres to (see the Transparent Methods for details).

Cells
A cell is described as a set of discrete lattice sites in a cellular Potts model (CPM) (see Figure 1A). Cells in the CPM change shape by iteratively making extensions and retractions of the cell boundary to adjacent lattices sites. This process also accounts for the formation and degradation of adhesions with the substrate. Time is measured in so-called Monte Carlo Steps (MCS), during which all sites along the cell boundary are given the opportunity to move. An additional rule is introduced to model the adhesive effect of FAs (assumption (6)): an energy barrier proportional to the size of the FA reduces the probability that the cell retracts from the substrate.

Cell Traction Forces
In addition to the retraction forces given by the CPM, we assume that the cell applies contractile forces onto each of the FAs. Based on the shape of the cell, we calculate the contractile force using the First Moment of Area (FMA) model (Lemmon and Romer, 2010) ( Figure 1B). This model assumes that the cell acts as a single contractile unit, resulting in a force onto each of the FAs pointing toward the center of mass and proportional to the distance to the center of mass. The forces at the individual FAs build up slowly from the latest force applied on the FA during the last MCS, to the force given by the FMA model at the given site and

OPEN ACCESS
iScience 23, 101488, September 25, 2020 3 iScience Article current MCS. For the latter, we use a model previously proposed by  in which the force builds up faster on stiffer substrates (assumption (5)).

Focal Adhesions
We model FAs as clusters of catch-slip bonds (assumption (2)), as proposed by . This model is based on experimental data of a single a5-b1 integrin that has constant binding rate (assumption (3)) and assumes that the degradation rate of bound integrins decreases with force (assumption (4)). The assembly and disassembly of a FA is then given by the collective dynamics of a cluster of integrin bonds. At each site of our 2D CPM, we implement the dynamics of one such FA (Figure 1C). Note that the FAs vanish near the cell center, because of insufficient traction. To increase visibility in the next figures, we visualize every other FA and show the average FA size of the surrounding four neighboring lattice sites.

Simulation
The simulations proceed as follows. In Model 1, we initiate cells in the CPM ( Figure 1A). Then, we let the forces build up for t FA seconds ( Figure 1B) and let the FAs grow simultaneously ( Figure 1C). After that, we let the cells move for one MCS ( Figure 1A) and repeat. Model 2 includes an extra step in which planar stress in the linear elastic ECM substrate (assumption (1)) reinforces the FAs (assumption (7)) (Model 2.1) or increases the forces applied to them (Model 2.2) ( Figure 1D).

Parameter Values
Where possible, parameter values were chosen based on literature values (Table S1), or else estimated, followed by parameter sensitivity studies.

Catch-Bond Cluster Dynamics Suffices to Predict Cell Area as a Function of Substrate Stiffness
We first tested whether the build-up of stiffness-dependent forces alone suffices to explain the observed cell deformation. We simulated Model 1 on a substrate of size 500 mm by 500 mm with stiffness of 1kPa, 5kPa, and 50kPa for 2000 MCS (z5.5 h) ( Figure 2A and Video S1). On the softest substrate of 1kPa, FAs did not grow and the cell did not spread. On a slightly stiffer substrate of 5 kPa, FAs grew and the cell increased in area. On the stiffest substrate of 50 kPa, the cell area further increased. Large FAs became visible around the cell perimeter ( Figure 2A). The cell area increased 2.5-fold (2500 mm 2 -6500 mm 2 ) as stiffness increased up to 50 kPa ( Figure 2B), consistent with experimental observations (Asano et al., 2017;Califano and Reinhart-King, 2010). We compared our model predictions with empirical functions relating substrate stiffness to cell area proposed in  and in  (Figures S1B-D) and found these to be consistent, apart from incidental quantitative discrepancies (See Table S2).
Comparing predicted and observed spreading kinetics ( Figure 2C), we found qualitative resemblance for endothelial cells Yeung et al., 2005) and fibroblasts Yeung et al., 2005). Cells spread to their final size over 30 min, a timescale consistent with fibroblasts (Yeung et al., 2005), but somewhat shorter than endothelial cells in vitro (100 min ). To compare the timescale of spreading, we consider t 50 : the time for a cell to reach 50% of its maximal area. To estimate t 50 , we fit an error function AðtÞ = A 50 , 1 + erf tÀt50 t  to the simulated curves AðtÞ. In our simulations, we found that t 50 z5 min (with 95% confidence interval of [5.02,5.46]) for cells on a 50kPa substrate ( Figure S1K). In , t 50 z50 min for BAECs on RGD-derivatized PA gel. The model parameter t FA sets the time scale of FA maturation (time for FA growth before the cell can extend or retract). Higher values of t FA result in slower cell spreading ( Figure S5A).
We next analyzed the dependence of final cell area on substrate stiffness for a range of values of t FA (Figure S5B). Because cell spreading and stiffness were most strongly correlated at t FA = 10 s (slightly less than the lifetime of cellular protrusions (Knorr et al., 2011)), we used this as a default value for other simulations.  For human fibroblasts , the timescale of spreading, t, decreases as the stiffness of the substrate increases ( Figure S1E). Fitting the inverse exponential function ( Figure S1J) from , we obtain a similar dependence of t on stiffness ( Figure S1H). As shown in Figure S1J, our fit is 4 times faster, but, as previously shown ( Figures S5G and S5H), varying free parameters can produce even closer agreement. In , the rate of change of the area ( dA dt ðtÞ) was shown to initially increase before dropping down to 0 ( Figure S1G). On the stiffest substrate, dA dt ðtÞ was overall much faster (Figure S1G). Interestingly, even without explicitly putting such dynamics into our model, a very similar curve emerges ( Figure S1I).
We next investigated the distribution of the FA sizes ( Figure 2D). At 1kPa, all FAs did not grow above the nascent adhesion threshold (N 0 = 50), so they were assumed to break apart (see Transparent Methods). The median cluster size and variance were unaffected by substrate stiffness, in contrast with Prager- Khoutorsky et al. (2011) and Trichet et al., 2012 who observed larger adhesions on stiffer substrates. Interestingly, and in line with our model observations,  report that average adhesion size is independent of substrate stiffness, although cell area and number of adhesions increase. A more detailed analysis of the distribution of the FAs revealed that stiffer substrates enhance the proportion of larger FAs, as observed in Paszek et al. (2005) (20% on 50,000 kPa, 15% on 10 kPa, and 10% on 5 kPa for FA with N>100). Furthermore, on stiffer substrates, large FAs were found not far from the cell periphery ( Figure 2A for 50 kPa and Figure S1). Similarly, in Shemesh et al. (2012) and in Shemesh et al. (2009) FAs were observed at the lamellipodium-lamellum interface, a small distance away from the cell edge. Choi et al. (Choi et al., 2008) attributed this observation to slow integrin-actin association, but in our model, this emerges from the fact that cell interior forces have more time to develop, thus stabilizing adhesions. On softer substrates, the largest FAs are found closer to the cell center (Figures 2A and S1). In the present model, we have assumed that the ECM can be approximated as a homogeneous, linearly elastic material. Although this can be a suitable approximation for small deformations of synthetic matrices, such as functionalized polyacrylamide matrices (Storm et al., 2005), natural matrices are highly inhomogeneous and nonlinearly elastic. To assess the effect of local inhomogeneities of ECM stiffness or of fluctuations in stiffness measurements by the cells (Beroz et al., 2017), we next simulated ECM inhomogeneity by applying random spatial variations to the Young's modulus of the ECM ( Figure S4A). Interestingly, the degree of inhomogeneity had little effect on the cell spreading area, because spreading is due to the average behavior of the FAs in the whole cell.
Mechanism of Uniform Cell Spreading. All in all, the mechanosensitive kinetics of the FAs suffices to predict stiffness-dependent cell area and spreading dynamics in a noise-resistant way. A schematic overview of our proposed mechanism for uniform cell spreading on soft versus stiff substrates is shown in Figure 2E. Cells (gray) send out protrusions (cyan) to probe the substrate by gradually building up forces. On a soft substrate, forces develop slowly, and thus adhesions do not grow enough to stick the cell to the substrate, and the cell retracts its protrusion. On a stiff substrate, forces build up rapidly, causing FAs to stabilize so the cell adheres to the substrate. This process repeats with new protrusions, allowing for maximal cell spreading.

Focal Adhesion Strengthening due to Matrix Stress Induces Cell Elongation
Model 1 correctly predicts the effect of substrate stiffness on cell spreading, but it cannot yet explain cell elongation. As cells pull, they deform the ECM, generating planar stresses. Uniaxial stretching of the substrate speeds up FA assembly, leading to FA strengthening and cell elongation, whereas on a longer timescale, planar stress induces FA disassembly, dampening cell elongation . We therefore hypothesize that an effect of planar stress on FAs may be involved in stiffnessinduced cell elongation.  iScience Article of the cell. On even more rigid substrates, the cells return to a circular shape. The same biphasic dependence of cell eccentricity on substrate stiffness was also experimentally observed for hMCS cells elongating most strongly on substrates of 10 kPa (Zemel et al., 2010). In general, the substrate rigidity associated with maximal elongation is cell-type and matrix-composition dependent. But since only a small range of substrate stiffness is usually tested, information is lacking about the exact stiffness at which cells start to elongate. For example, fibroblasts do so at 2kPa on collagen-coated PA gels. On fibronectin-coated PA gels, fibroblasts with PTK knockdown failed to elongate strongly at 30kPa but did so at 150kPa (Prager- Khoutorsky et al., 2011). MCS elongated at 9kPa (but not at 0.7kPa (Rowlands et al., 2008)), endothelial cells at 1kPa (Califano and Reinhart-King, 2010), and cardiomyocytes at 5kPa (Chopra et al., 2011).
We next tested how local inhomogeneities in substrate stiffness affect the ability of cells to elongate. As before, we added uniform noise to the substrate stiffness. Figure S4B shows the cell eccentricity as a function of ECM stiffness and the level of inhomogeneity; together these data suggest that the observed cell elongation is robust to noise. Next, we investigated the effect of inhomogeneities on a longer length scale, by imposing a sinusoidal function of substrate stiffness and varying its period and amplitude (see Transparent Methods and Figures S4D-L). Cells elongate as before, unless the period of the sinusoid is much greater than the cell diameter and the stiffness significantly deviates from its baseline (50kPa, Figures  S4K -L). This suggests that if the ECM stiffness changes significantly on a length scale of a cell, the cell cannot find two stiff-enough anchor spots so as to elongate.
We compared the distribution of FAs of elongated cells to uniform cell spreading. Figure 3C shows the distribution of the FA sizes as a function of substrate stiffness. As before, the median cluster size does not vary much with substrate stiffness. The shape of the distributions, however, are flatter, with higher variance for elongated cells compared with round cells (kurtosis kz2:0, SD z2700 for elongated cells and kz 2:3, SD z 1500 for round cells). This larger variation in FAs sizes in elongated cells stems from the polarized force field that forms large FAs at the poles and small FAs at the lateral sides.
We next performed a sensitivity analysis for parameters not constrained by experimental data. For increased values of p, which regulates planar stress-induced FA strengthening, cells elongate over a much wider range of stiffnesses ( Figures S2A-B Because of the relatively large cellular temperature T, the cell contours are fairly rugged. This temperature value was required to stabilize the elongated cell shapes. Although a lower temperature can still produce similar behavior ( Figure S6A), we observed that elongated cells go through cycles of collapse (when contractile force is higher than adhesive force) and elongation. Higher temperature ensures that cells are sufficiently motile, rapidly forming new protrusions and adhesions after a slight retraction.
Model 2.2: Planar Stress Enhances Traction Force. We also tested an alternative model in which the force exerted on the FAs increases in response to planar stress. We assumed that the stall force increases as a function of matrix stress, i.e., . This alternative mechanism produces ECM-stiffness-dependent cell morphologies that are similar to the default model ( Figure S3 and Video S3). This alternative mechanism could have various molecular origins. For instance, addition of vinculin through talin stretching can induce increased traction forces (Dumbauld et al., 2013). Stretching forces also induce a-smooth muscle actin recruitment to stress fibers (Goffin et al., 2006) and myosin motor binding (Uyeda et al., 2011).

Mechanism of Cell Elongation.
In conclusion, our model suggests a mechanism for cell elongation that we illustrate with a series of snapshots in Figure 3D. Because of stochastic variations, slightly eccentric cell shapes are formed occasionally. The first-moment-of-area (FMA) model (Lemmon and Romer, 2010)  iScience Article circles in Figure 3D1,2) in the underlying substrate. Because we assumed that ECM stress stabilizes FAs (Model 2.1), protrusions around such hot-spots stick to the ECM. Then, forces can further develop, ECM stress increases, and FAs further stabilize. As the cell elongates, the force field also polarizes ( Figure 3D3). On the lateral sides of the cell, traction forces are smaller and thus insufficient to maintain stable FAs (Figure 3D4). Together with the positive feedback between stress and FAs, the retractions at the lateral sides reinforce cell elongation on substrates of intermediate stiffness. In version 2.2 of our model, cell elongation happens similarly, but here ECM stress positively feeds back on cell traction forces, rather than on the adhesive force. On soft substrates, force build-up is not strong enough to stabilize the FAs, whereas on stiff substrates, all FAs stabilize, leading to round cells on both soft and stiff substrates.
Contractility Rate Changes Stiffness Regime on which Cells Elongate. There are large differences in the response to ECM stiffness between cell types. Fibroblasts spread out further on stiff ECMs than on soft ECMs, whereas endothelial cells do so to a lesser extent (Yeung et al., 2005). Cells display a large variability in contractility (Murrell et al., 2015) and in the actomyosin contraction rate, which may be suppressed by microtubules (Bershadsky et al., 2003a;Dugina et al., 2016). To test whether variability in the response to ECM stiffness could be due to differences in the rate of actomyosin fiber contraction, we varied cell contractility by modulating the velocity of the myosin motors, v 0 , during force build-up (see Equation 5 in Transparent Methods). We studied a range from 10 nm/s, corresponding to nonmuscle myosin-IIB  to 1,000 nm/s, corresponding to muscle myosin (Vogel et al., 2013). Figures 4A and 4B show the cell configurations for reduced (v 0 = 10 nm=s) and increased motor velocities (v 0 = 1000 nm=s), corresponding to velocities one order of magnitude below and above the default value (v 0 = 100 nm=s). Figures  4C and 4D plot the cell area and cell eccentricity as a function of motor protein velocity averaged over 25 independent simulations. With reduced motor velocity, the cells spread only on the stiffest matrices tested ( Figure 4C) and failed to elongate ( Figure 4D). Cells with increased motor velocities, however, elongated at relatively soft matrices of 5 kPa ( Figures 4B and 4D). We next asked whether our model mechanism also suffices to predict durotaxis: cell migration up a stiffness gradient (Lo et al., 2000). It has previously been proposed that cells durotact due to increased FA maturation on stiff parts of the matrix, whereas the FAs disintegrate and cellular protrusion retract at softer parts of the substrate (Wong et al., 2014). We tested this in-silico by placing cells on stiffness gradients and following their movement for 10,000 MCS ( z28hÞ. The stiffness increased linearly over the x axis from 1 kPa (left) to 26 kPa (right), corresponding to a slope of 20 Pa/mm. The cell starts at ð250 mm; 250 mmÞ at 6 kPa. We found that cells gradually drifted toward the right, the stiffer part of the gradient on homogeneous matrices (Video S4).
The simulated cells move in the x-direction at constant velocity ( Figure 5B), with an average speed of 4:3 mm=h, (average slope, dx=dt, over 25 simulations), comparable to 6.2 mm/h for mesenchymal stem cells Mechanism of Durotaxis. We illustrate our proposed mechanism for durotaxis in Figure 5F. The cell probes the substrate by randomly sending out protrusions. Because forces build up faster on the stiff side of the substrate, FAs grow larger than on the soft side. Protrusions are thus more likely to stick to the stiff part, allowing the cell to move forward as it retracts its back and new protrusions successfully stick in front. We tested the role of the feedback between ECM stress and FA stabilization, by setting p = 0, equivalent to running Model 1. Interestingly, cells durotact (but slightly slower, at 3:9 mm=h), suggesting that the catch-bond integrin dynamics are sufficient for durotaxis and any further stabilization of FA proteins makes cells move quicker.
The durotaxis mechanism requires that cells are sufficiently flexible to make new protrusions and retract old ones. Indeed, increasing the value of l that decreases cell flexibility (See Equation 3  Durotaxis speed should also depend on how reliably matrix stiffness controls the rate of protrusions and retractions. In the CPM, the stochasticity of such shape changes is controlled by the ''cellular temperature'' T (See Equation 2 in Transparent Methods), where for larger values of T, retractions are less tightly controlled by cell-substrate adhesion and other environmental and cellular properties. As predicted, cell speed is reduced for larger values of the cell motility parameter ( Figure 5C). We also tested to what extent matrix inhomogeneity affects the durotactic migration speed of cells; even noise levels of G 10 kPA did not affect the durotactic cell velocity ( Figure S4C).

DISCUSSION
We have shown that mechanoregulation of the assembly and disassembly of FAs by matrix stiffness and planar stress suffices to explain (a) increased cell spreading on stiff substrates (Figures 2A-B), (b) cell elongation on substrates of intermediate stiffness ( Figures 3A-B), and (c) durotaxis ( Figure 5A). Many previous mathematical models of cell spreading and cell migration used a bulk equation for cell-matrix adhesion, motivated by the underlying FA dynamics (Bischofs and Schwarz, 2003;Kabaso et al., 2011;Ni and Chiang, 2007;Ramos et al., 2018;Zemel et al., 2010). More closely related to the present work, hybrid cell and FA models were developed to study how the mechanosensitive growth of FAs direct cell spreading and migration (see, e.g., Copos et al. (2017) and Uatay (2018)). It was proposed that cell spreading is regulated by a matrix-stiffness-induced stress fiber persistence and upregulated traction forces (Ronan et al., 2014;Vernerey and Farsad, 2014). Another model has proposed an enhanced cell contraction on stiff matrices (Stolarska and Rammohan, 2017). Because this effect counteracted cell spreading, the model could not explain increased cell spreading on stiff matrices. These previous models suggest that an upregulation of cell forces is required to induce cell spreading on stiff matrices. Similarly, Novikova and Storm  noted that force evolution alone does not enable a cell to distinguish soft from stiff matrices, because ultimately the same stall force is reached. They propose that stiffness sensing requires cells to apply more force on stiffer substrates. Our model suggests that force evolution can be sufficient for cells to sense matrix stiffness. In our model, cells can more rapidly build up forces ll OPEN ACCESS iScience 23, 101488, September 25, 2020 11 iScience Article on stiffer matrices during the lifetime of a cell protrusion, allowing FAs to stabilize and protrusions to stick to the matrix and subsequently the cell to spread. On soft matrices, the force does not reach a sufficient level during the lifetime of a protrusion to allow for stabilization. So, force evolution together with the dynamic nature of protrusions/retractions enable the cell to distinguish stiff from soft.
For uniform cell spreading, it is essential that (a) FAs apply an adhesive force on the cell (clearly the case, based on the known function of FAs), and (b) the rate of force build-up depends on ECM stiffness. Otherwise, cells are not able to distinguish between soft and stiff matrices and spread accordingly. After adding the assumption that planar stress contributes to the stabilization of FAs, the predictions of our model also sufficed for cell elongation. A possible molecular mechanism for this effect is the stretching of talin in FAs.
Stretching of talin exposes vinculin-binding sites (Rio et al., 2009). Vinculin in turn binds the FA to the cytoskeleton, which strengthens cell-matrix adhesion . In agreement with this hypothesis, vinculin regulates cell elongation on glass substrates (Ezzel et al., 1997). Although observations in fibroblasts suggest that traction forces are independent of ECM stiffness (Feld et al., 2020;Freyman et al., 2002), as we assumed in this model, there are also contradicting observations: (1) vinculin increases cell traction forces (Dumbauld et al., 2013), and (2) stressing FAs induces a-smooth muscle actin recruitment to stress fibers, which in turn increases traction forces (Goffin et al., 2006). We have, therefore, also tested an alternative mechanism in which planar stress induces an increase in cell traction forces ( Figure S3); this mechanism also suffices for cell elongation. Future, quantitative comparisons of our model predictions and in vitro observations may help elucidate which of these two mechanisms, if any, best explains cell elongation and for what cell types.
Our model also predicts distinct ranges of substrate stiffness for cell elongation could be due to diverse myosin motor velocities (Figure 4), possibly accounting for cell-type dependence (Georges and Janmey, 2005;Yeung et al., 2005). For example, in ovarian tumor cells deficient in Dlc1 (responsible for phosphorylation of nonmuscle IIA myosin (Sabbir et al., 2016)), elongation is promoted, suggesting that an increase in motor velocity indeed facilitates greater cell elongation. This prediction is experimentally testable by studying cells that express different isoforms of myosin or by overexpression or inhibition of a given myosin isoform in cells while systematically varying substrate stiffness.
The assumptions of our model also sufficed for durotaxis. The mechanoresponsivity makes FAs longerlived on the stiffer side of the matrix than on the softer side of the matrix. The resulting bias in FA turnover and pseudopod turnover is responsible for a drift of the cell toward the stiffer side of the matrix. Instead of such implicit, emergent effects, previous models of durotaxis have often assumed direct effects of ECM stiffness on FA density and polarity (Kim et al., 2018;Yu et al., 2017). Feng et al. (Feng et al., 2018) showed that if FA degradation is higher in the back than in the front and FAs mature under applied force, then a cell can durotact. Based on experimental observations, Novikova et al. assumed that cells move more persistently on stiffer substrates and showed that a persistent random walk can reproduce durotaxis (Novikova et al., 2017). In contrast to the direct effect of ECM stiffness on cell polarization, in our model durotaxis emerges from a bias in FA turnover. By knocking-out the model feedback between ECM stress and FA stabilization (obtaining Model 1), we found that durotaxis is still possible, but at a slightly reduced speed. This suggests that the catch-bond integrin behavior suffices for durotaxis. The role of vinculin could then be to increase cell speed. This prediction can be tested empirically by knocking down vinculin or related structural proteins.
Our model predicts that durotaxis speed saturates with steeper stiffness gradients. In our model, this emerges from the fact that the growth of FAs saturates at high ECM stiffness, giving an upper bound to the cell speed. To validate this prediction, a wider range of stiffness gradients should be investigated.

Model Limitations
No model can account for every experimental observation and ours is no exception. Some experimental data are clearly inconsistent with our predictions. For example, human cancer cells move faster on stiffer substrates, whereas our simulated cells are slower. For smooth muscle cells, it was reported that cell velocity has a biphasic dependence on substrate stiffness, i.e. cells move slowest on the softest and stiffest substrates (Peyton and Putnam, 2005 iScience Article biochemical signaling acts to shape, regulate, and fine-tune the cell's mechanical sensing and response, aspects that are not yet considered in our model. A bridge between the mechanical and biochemical signals for cell migration is a promising avenue for future models. For computational simplicity, our model represents the ECM as a uniform, isotropic, linearly elastic material (e.g., synthetic polyacrylamide gel), with fiber sizes that are sufficiently small relative to the size of the cell. So, our predictions are most accurate for cells moving on PA gels. In natural matrices, inhomogeneities result from traction-force-induced fiber realignment, density changes, and nonlinear strain-stiffening. As a first step toward more complex matrices, we have studied how cell spreading, elongation, and durotaxis are affected by inhomogeneity in substrate stiffness (uniform noise, Figures S4A-C) and explored the effect of small-scale spatial inhomogeneities on cell elongation ( Figures S4D-L). In our ongoing work, we are extending the model to better represent the fibrous and strain-stiffening properties of natural ECMs (Hall et al., 2016;Han et al., 2018;van Helvert and Friedl, 2016) by incorporating discrete ECM models (Feng et al., 2015(Feng et al., , 2018Kim et al., 2017;Licup et al., 2015) with our cell-based models. Although the effect of traction forces propagates further into nonlinearly elastic, fibrous ECMs (Ma et al., 2013), the key results of our model are based on increased stability of FAs on strained matrices and so, will likely still apply.
Our model currently does not accurately predict that FAs are larger on stiff than on soft ECM substrates ( Figure 2D) (Prager-Khoutorsky et al., 2011). Although there is a small increase in size in the range between 1 kPa and 10 kPa, for large stiffnesses the FAs do not increase in size. This might be due to the assumed fixed pool of free integrin bonds that reduces the growth rate of new FAs once many FAs already exist. Furthermore, our lattice-based model does not consider spatial, cooperative effects in integrin clustering: small clusters may merge into larger adhesions. Integrins diffuse along the cell membrane and they are activated upon interaction with regulatory proteins such as talin (Welf et al., 2012) and vinculin (Humphries et al., 2007). Similar mechanisms include phosphorylation of p130cas in response to stretching, which activates the small GTPase Rap1 (Sawada et al., 2006) and resulting integrin activation (Bos et al., 2003). Future versions of our model will consider these and further regulatory mechanisms of FA dynamics (Ali et al., 2011;Deshpande et al., 2008;Vernerey and Farsad, 2014;Welf et al., 2012) and thus refine its predictive value.
Another simplifying assumption concerns the distribution of the cellular traction forces (Lemmon and Romer, 2010). Although this model has been experimentally validated, it has some limitations, principally, that it is inconsistent with the Hamiltonian-based CPM forces. (The Hamiltonian acts as a potential energy, whose gradient defines a force field, at least along the cell edge; see Rens and Edelstein-Keshet (2019).) Alternative force descriptions have also been proposed in . These formulations also lead to different spatial distributions of forces inside the cell. In Rens and Edelstein-Keshet (2019), for example, it was shown that linear, polynomial, or exponential dependence of force on distance from the cell centroid were all consistent with experimental traction force data. In real cells, traction force distributions depend on the details of the stress fiber localizations and substrate topography (Soiné et al., 2015). In models that we have recently developed (Pomp et al., 2018;Schakenraad et al., 2020), cell geometry affects cytoskeletal orientation, which in turn determines the directionality of stress fibers. Integrating these new traction force models into the CPM will result in future improvements and greater accuracy.

Limitations of Study
Our work is computational, and, although based on experimental literature, has yet to be (in)validated against dedicated future experimental studies. The model assumes a homogeneous linearly elastic ECM, among many simplifications. Hence, it falls short of explaining biological behavior that stems from the true fibrous, nonlinear elastic nature of the ECM. The model presently omits intracellular signaling cascades that sense and respond to mechanochemical stimuli and thus, in its current form, does not accurately account for every experimentally observed cell behavior nor is it calibrated to specific cell types. To apply our model to more specific cell systems, additional cellular and adhesion mechanisms and cell-extracellular matrix interactions should be considered. Finally, any modeling platform has limitations as well as benefits. Although the cellular Potts framework allows for highly resolved cell shape computations, it comes with drawbacks, and does not explicitly describe cell forces nor inherent cellular dynamics time scales.

Materials Availability
There are no materials associated with this work.

Data and Code Availability
The model (C++ code) and parameter files discussed in the paper are available from GitHub at https:// github.com/rmerks/FA-CPM-FEM. The code was written in the Cellular Potts-Finite Element framework described in .

METHODS
All methods can be found in the accompanying Transparent Methods supplemental file.
van Helvert, S., Storm, C., and Friedl, P. (2018).      Model results for lower temperature (T = 0.5). The spreading strength parameter λ C was increased to 1000, otherwise cells would not spread. Notice that at this lower temperature, an elongated cell shape is unstable (see the three panels to the right): the cell repeatedly rounds up and elongates again. (B) Rescaling the grid by a factor of s = 2: ∆x = ∆x/s, λ = λ/s 2 , J = J/s, (C) Rescaling the grid by a factor of s = 4. In (C), cells disappear at 1kPa (they become 1 pixel in size). In both (B) and (C), the elongated shape at 50kPa is unstable like in (A).

Transparent Methods
We developed a multiscale model where cell movement depends on force induced focal adhesion dynamics. The model couples a cell-based model, substrate model and focal adhesion model in the following way. The Cellular Potts Model (CPM) describes cell movement. The shape of the cell is used to describe the stall forces that the cell exerts on the focal adhesions attached to a flexible substrate. These forces affect the growth of the focal adhesions. We assume that focal adhesions are clusters of integrins that behave as catch-slip bonds. Its dynamics are described using ordinary differential equations (ODEs). Finally, we assume that the cell-matrix link is strengthened by matrix stresses, which we calculate using a finite element model (FEM). The default parameter values are described in Table S1.

Cellular Potts Model
To simulate cell movement, we used the Cellular Potts Model (CPM) (Graner and Glazier, 1992). The CPM describes cells on a lattice Λ ⊂ Z 2 as a set of connected lattice sites. Since the simulations in this article are limited to one cell, we describe the CPM here for a single cell. To each lattice site x ∈ Λ a spin s( x) ∈ {0, 1} is assigned. This spin value indicates if x is covered by the cell, s( x) = 1, or not, s( x) = 0. Thus the cell is given by the set, (1) The cell set C evolves by dynamic Monte Carlo simulation. During one Monte Carlo Step (MCS), the algorithm attempts copy a spin value s( x) from a source site x into a neighboring target site x from the usual Moore neighbourhood. Such copies mimic cellular protrusions and retractions. During an MCS, N copy attempts are made, with N the number of lattice sites in the lattice. Whether a copy is accepted or not depends on a balance of forces, which are represented in a Hamiltonian H. A copy is accepted if H decreases, or with a Boltzmann probability otherwise, to allow for stochasticity of cell movements: (2) Here ∆H = H after − H before is the change in H due to copying, and the cellular temperature T ≥ 0 determines the extent of random cell motility. Furthermore, Y denotes a yield energy, an energy a cell needs to overcome to make a movement. Finally, to prevent cells from splitting up into disconnected patches, we use a connectivity constraint that always rejects a copy if it would break apart a cell in two or more pieces. Following , we use the following Hamiltonian: (3) The first term of H denotes cell contraction, where A is the area of the cell and λ is the corresponding Lagrange multiplier. In the second term, J(s( x), s( x )) is the adhesive energy between two sites x and x with spins s( x) and s( x ). When taking a sufficient large neighborhood, this term describes a line tension, as it approximates the perimeter of a cell (Magno et al., 2015). We take a neighborhood radius of 10 for this calculation. The third term describes the formation of adhesive contacts of cells with the substrate, where the bond energies lower the total energy , causing the cells to spread. The parameter λ C is the corresponding Lagrange multiplier. The energy gain of occupying more lattice sites saturates with the cell area, as the total number of binding sites is limited. The parameter A h regulates this saturation.
To describe cell-matrix binding via focal adhesions, we implement the following yield energy in the CPM where N ( x ) is the size of the focal adhesion at the target site. This models that a retraction is energetically costly for a cell to make, because it needs to break the actin-integrin connection. We assume that the size of the actinintegrin link is proportional to the size of the focal adhesion, i.e. the number of integrin bonds (Boettiger, 2007), and that the strength of a focal adhesion saturates  with a parameter N h . The substraction of N 0 represents that a focal adhesion only creates extra linkage if it is greater than a nascent adhesion. Note that the Y can not become negative, because we assume that focal adhesions smaller than N 0 , a nascent adhesion, breaks down due to its short lifetime, see section 'Focal Adhesions'. So, only focal adhesions larger than N 0 create a yield energy. In section 'Substrate stresses', we further adapt this yield energy to describe a matrix stress induced focal adhesion reinforcement.

Cell traction forces
Following , we assume that traction forces are generated by myosin molecular motors on the actin fibers, of which the velocity is given by where v 0 is a free velocity. The traction forces are applied to the ECM, which we assume is in plane stress. The constitutive equation is given by h ∇σ = F where σ is the ECM stress tensor and h is the thickness of the ECM. We assume that the ECM is isotropic, uniform, linearly elastic and we assume infinitesimal strain theory. We solve this equation using a Finite Element Model (FEM) (see section 'Substrate stresses'). In the FEM, traction field f and ECM deformation u are related by: where K is the global stiffness matrix given by assemblying the local stiffness matrices K e for each lattice site e where B is the conventional strain-displacement matrix for a four-noded quadrilateral element (Davies, 2011) and E is the Young's modulus and ν is the Poisson's ratio of the ECM. For more details on this part of the model, we refer to our previous work (van . Following , the force build-up is given by the ODE: The matrix K describes force interactions between neighbouring nodes in the ECM. However, since solving this equation is expensive, we ignore the interactions between neighbouring sites, i.e., we reduce K to a scalar for each site x. This gives us, where F 0 is the force already exerted by the actin and t k = | Fs| v0K . Here, K is given by the diagonal entry of K at site x, i.e. the stiffness of this node, neglecting changes in local stiffness due to connections to neighbouring nodes described in the off-diagonal entries of K. Since the cell configuration and therefore the traction forces change each MCS, the tension on the focal adhesions does not build up from zero, but from the tension that was built up during the previous MCS: F 0 at the current MCS is given by F (t FA ) of the previous MCS.
To calculate the stall force of the actin fibers, F s , we employ the empirical first-moment-of-area (FMA) model (Lemmon and Romer, 2010). This model infers the stall forces from the shape of the cell of the CPM, based on the assumption that a network of actin fibers in the cell acts as a single, cohesive unit, So, the force at site x is calculated as the sum of forces between x and all other sites y within cell C that are connected to x (this sum excludes line segments [ y x] running outside the cell that occur if the shape of cell C is non-convex).
The force is assumed to be proportional to the distance between the sites. We divide over the cell area A such that force increases roughly linear with cell area, as experimentally observed (Califano and Reinhart-King, 2010b).

Focal adhesions
At each lattice site occupied by the cell, x ∈ C, a focal adhesion is modeled as a cluster of bound integrin bonds N . Each individual integrin bond behaves as a catch-slip bond, whose lifetime is maximal under a positive force . Accordingly, the growth of a cluster of such bonds is described by the ODE-model derived by , with γ is the binding rate of integrins to the ECM, N a the number of free bonds, and N b the maximum number of bound bonds a lattice site can contain. Following the previous work , we assume that the closing of integrin bonds occurs at a constant rate once the first adhesion bonds bring the cell membrane in sufficiently close proximity to the ECM; we assume that this is the case when the cell covers the lattice site. The logistic growth term, γN a (t) 1 − N ( x,t) , is a slight modification of the previous work  and avoids packing more integrins in a lattice site than the size of the lattice site (∆x 2 ) can accommodate. The degradation of the focal adhesions d(φ) depends on the tension φ on the focal adhesion N . This degradation rate is given by where φ s and φ c are nondimensional parameters that describe the slip and catch bond regime, respectively. Parameter d 0 is the baseline unbinding rate of the integrins . Here, φ( x, t) = | F ( x,t)| ∆x 2 is the stress applied to the lattice site of the focal adhesion. To conform the units of force φ( x, t) to the dimensionless parameters φ s and φ c , we scale it with α (units #integrins N/m 2 ) in Equation 12. We assume that the number of free bonds N a is limited by the number of available integrin receptors in the entire cell, N m . These N m receptors can be recruited to each focal adhesion site and enable binding of a bond. Thus, N a (t) = N m − x∈C N ( x, t). We let the focal adhesions grow after each MCS for t FA seconds with time increments of ∆t FA . If there is no pre-existing focal adhesion at site x ∈ C, we set N ( x) = N 0 , so that at this site, a new initial adhesion is formed. This assumption represents the generation of focal complexes or nascent adhesions, precursors of focal adhesions that contain a small amount of integrins and have a very short lifetime (Changede and Sheetz, 2016). Also, after the focal adhesions were allowed to grow, i.e. after t = t FA seconds, we set all N ( x) < N 0 back to N ( x) = N 0 , again modeling the quick (re)generation of nascent adhesions.
When a site x is removed from the cell C after a retraction, such that s( x) = 0, we set N ( x) = 0 reflecting the destruction of the focal adhesion. We assume that if a cell extends, i.e. a site x is added to the cell C, a nascent adhesion is formed: we set N ( x) = N 0 .

Substrate stresses
The forces that were build up during a MCS, F (t FA ) are applied as planar forces to a finite element model (FEM). The FEM calculates the stress tensor σ( x) on each lattice site. We assume that the integrin-cytoskeletal adhesion strengthens as a result of stress. We define g(σ) = 1 2 (σ xx + σ yy ) if 1 2 (σ xx + σ yy ) ≥ 0 0 if 1 2 (σ xx + σ yy ) < 0 the positive hydrostatic stress of the stress tensor that describes how much stress the focal adhesion experiences.

Model 2.1
We extend the yield energy as follows: We thus assume that stress strengthens the focal adhesions, with parameter p and that this strengthening saturates with parameter σ h .

Model 2.2
In model 2.1, we proposed that matrix stress induces focal adhesion strengthening but noted that matrix stress might also reinforce cell contractility. Supplementary Figure3 shows the results of having F s = F s · 1 + p g(σ( x )) σ h +g(σ( x )) instead of focal adhesion strengthening as described in the main text. Since matrix stresses are defined on the lattice sites while forces are defined on the nodes of the lattice, we needed to assume some interpolation. We choose to take F s = F s · 1 4 surrounding4nodes 1 + p g(σ( x )) σ h + g(σ( x )) .

Stiffness gradient
To study durotaxis, we model a stiffness gradient in the x-direction on a lattice of 1250 µm by 500 µm. The Young's modulus of the substrate E(P a) is given by E(x) = max{1, 6000 + (x − 250) · slope}, with x in µm, such that the Young's modulus at the center of the cell at time t = 0 is 6000 Pa and is nonzero. The default value for the slope is 20 Pa/µm.

Noise in substrate stiffness
To simulate inhomogeneity in substrate stiffness, we add random noise between −E noise and +E noise to the Young's modulus of the substrate E at each lattice site. If this sets the substrate stiffness to a negative value, we cut it off at zero. We tested the effect of noise on cell spreading, cell elongation and durotaxis for values of E noise of 100 Pa, 1 kPa, 2 kPa, 5 kPa, 7.5 kPa and 10 kPa. The results are in Supplementary Figure4.

Local variations in substrate stiffness
To simulate longer range inhomogeneities in substrate stiffness, we modify the substrate stiffness as follows: where E A is the amplitude of the sinusoidal and T E its period. The results are in Supplementary Figure4.

Grid resolution
The spatial resolution of the cell and the number of integrins per CPM pixel can be changed by scaling ∆x. In Figure S6B and C, we show cell elongation for rescaling by a factor of 2 and 4 respectively. The caption indicates the rescaling of the parameters. The results are qualitatively similar, although cells start to elongate at softer substrates and at 50kPa, the elongated shape is unstable (as in Figure S6A). On 1kPa, the cell fails to spread at all for a grid scaling factor of 2 ( Figure S6B). Finally, the pool of integrins is also more quickly depleted, resulting in less focal adhesions near the cell edge. This causes instability of the elongated shape. When a cell is very long, its inward contractile force cannot be counteracted by the focal adhesion far from the cell edge.