Mechanical cell-matrix feedback explains pairwise and collective endothelial cell behavior in vitro

In vitro cultures of endothelial cells are a widely used model system of the collective behavior of endothelial cells during vasculogenesis and angiogenesis. When seeded in an extracellular matrix, endothelial cells can form blood vessel-like structures, including vascular networks and sprouts. Endothelial morphogenesis depends on a large number of chemical and mechanical factors, including the compliancy of the extracellular matrix, the available growth factors, the adhesion of cells to the extracellular matrix, cell-cell signaling, etc. Although various computational models have been proposed to explain the role of each of these biochemical and biomechanical effects, the understanding of the mechanisms underlying in vitro angiogenesis is still incomplete. Most explanations focus on predicting the whole vascular network or sprout from the underlying cell behavior, and do not check if the same model also correctly captures the intermediate scale: the pairwise cell-cell interactions or single cell responses to ECM mechanics. Here we show, using a hybrid cellular Potts and finite element computational model, that a single set of biologically plausible rules describing (a) the contractile forces that endothelial cells exert on the ECM, (b) the resulting strains in the extracellular matrix, and (c) the cellular response to the strains, suffices for reproducing the behavior of individual endothelial cells and the interactions of endothelial cell pairs in compliant matrices. With the same set of rules, the model also reproduces network formation from scattered cells, and sprouting from endothelial spheroids. Combining the present mechanical model with aspects of previously proposed mechanical and chemical models may lead to a more complete understanding of in vitro angiogenesis.

naling during morphogenesis, and introduce a new set of computational tools for modeling mechanical interactions between cells and the extracellular matrix. inhibition of chemotaxis [24] are needed to transform these circular aggregates into vascular-like network patterns. Related network formation models studied the role of ECM-bound growth factors [25][26][27] and a range of additional secreted and exogenous growth factors [27], and studied the ability of the contactinhibition mechanism to produce three-dimensional blood-vessel-like structures [28]. Szabó and coworkers found that in culture, astroglia-related rat C6 cells and muscle-related mouse C2C12 cells organize into network-like structures on rigid culture substrates [29], such that ECM-density or chemoattractant gradients are excluded. They proposed a model where cells were preferentially attracted to or preferentially adhered to locally elongated structures. As an alternative mechanism for "gel-free" network formation it was found that elongated cells can also produce networks in absence of chemoattractant gradients [30].
Paradoxically, despite the diverse assumptions underlying the mathematical models proposed for vascular network formation, many are at least partly supported by experimental evidence. This suggests that a combination of chemotaxis, and chemical and mechanical cell-ECM interactions drives network formation, or that each alternative mechanism operates in a different tissue, developmental stage, or culture condition. A problem is that one mathematical representation may represent a range of equivalent alternative underlying mechanisms. For example, a model representing cell-cell attraction cannot distinguish between chemotaxis-based cellular attraction [20,21,23,24], attraction via haptotaxis [19], direct mechanical attraction [15,31] or cell shape dependent adhesion [29,32], because the basic principles underlying these models are equivalent [12,24]. As a solution to this problem, a sufficiently correct complete description of endothelial cell behavior should suffice for the emergence of the subsequent levels of organization of the system, an approach that requires that the system has been experimentally characterized at all levels of organization.
The role of cell traction and ECM mechanics during in vitro angiogenesis have been characterized experimentally particularly well, making it a good starting point for such a multiscale approach. Endothelial cells apply traction forces on the extracellular matrix, as demonstrated by a variety of techniques, e.g., wrinkle formation on elastic substrates [9], force-generation on micropillar substrates [33], and traction force microscopy [6,34]. Using scanning electron microscopy, Vernon and Sage [9] found that ECM ribbons radiate from endothelial cells cultured in Matrigel, suggesting that the traction forces locally reorient the extracellular matrix. The cellular traction forces produce local strains in the matrix, which can affect the motility of nearby cells [2]. Thus endothelial cells can both generate, and respond to local strains in the extracellular matrix, suggesting a feedback loop that may act as a means for mechanical cell-cell communication [2] and hence coordinate collective cell behavior. Here, we use a hybrid cellular Potts and finite element model to show that a set of assumptions mimicking mechanical cell-cell communication via the ECM suffices to reproduce observed single cell behavior [35,36], pairwise cell interactions [2], and collective cell behavior: network formation and sprouting.

Response of endothelial cells to static strains in ECM
First we set out to capture, at a phenomenological level, the response of endothelial cells to static strains in the ECM in absence of cellular traction forces. When grown on statically, uniaxially stretched collagen-enriched scaffolds, murine embryonic heart endothelial (H5V) cells orient in the direction of strain, whereas cells grown on unstrained scaffolds orient in random directions [37]. Because the collagen fibers make the scaffold stiffen in the direction of strain, we hypothesized that the observed alignment of cells is due to durotaxis, the propensity of cells to migrate up gradients of substrate rigidity [38] and to spread on stiff substrates [39,40]. In our model we assumed (a) strain stiffening: a strained ECM is stiffer along the strain orientation than perpendicular to it, such that (b) due to durotaxis the endothelial cells preferentially extend pseudopods along the strain orientation, along which the ECM is stiffest, giving cells the most grip. To keep the ECM mechanics simulations computationally tractable, we assumed an isotropic and linearly elastic ECM. With these assumptions it is not possible to model strain stiffening explicitly. We therefore mimicked strain stiffening by assuming that stiffness is an increasing, linear function of the local strain.
Durotaxis was modelled as follows, to reflect the observation that focal adhesion maturation occurs under the influence of local tension [41]: At low local stiffness, we applied standard cellular Potts dynamics to mimic the iterative formation and breakdown of ECM adhesions, producing "fluctuating" pseudopods. However, if the stiffness was enhanced locally, we assumed that the resulting tension in the pseudopod led to maturation of the focal adhesion [41,42], stabilizing the pseudopod as long as the tension persists. To mimic such focal adhesion maturation in the cellular Potts model, we increased the probability of extension along the local strain orientation, and reduced the probability of retraction (see Methods for detail). Figure 1 A shows the response of the simulated cells to uniaxial stretch along the vertical axis. With increasing values of the durotaxis parameter λ durotaxis (see Eq. 8), the endothelial cells elongate more. To test the sensitivity of the durotaxis model for lattice effects, we varied the orientation of the applied strain over a range [0−180] • and measured the resulting orientation of the cells. Figure 1 shows that the average orientation of the cells follows the orientation of the stretch isotropically. Thus the durotaxis component of our model phenomenologically reproduces published responses of endothelial cells to uniaxial stretch [37].

Generation of strains in ECM due to cellular traction
We next attempted to mimic the forces applied by cells onto the extracellular matrix, in absence of durotaxis. Traction-force microscopy experiments [34,39] show that endothelial cells contract and exert tensional forces on the ECM. The forces are typically directed inward, towards the center of the cell, and forces concentrate at the tips of pseudopods. A recent modeling study by Lemmon and Romer [43] found that an accurate prediction of the direction and relative magnitudes of these traction forces within the cell can be obtained by assuming that each lattice node i covered by the cell pulls on every other node the cell covers, j, with a force proportional to their distance, d i,j . Because this model gives experimentally plausible predictions for fibroblasts, endothelial cells, and keratocytes [43], we adopted it to mimic the cell-shape dependent contractile forces that endothelial cells exert onto the ECM. Figure 2 shows the contractile forces (black) and resulting ECM strains (blue) generated in our model by two adjacent cells. The traction forces and ECM strains become largest at the cellular "pseudopods", qualitatively agreeing with traction force fields reported for endothelial cells [34].
Mechanical cell-ECM feedback qualitatively reproduces effect of substrate stiffness on cell shape and motility The two previous sections discussed how the simulated cells can respond to and induce strain in the ECM in an experimentally plausible way. To test how the simulated cells respond to the strains they generate themselves, we studied the behavior of simulated, single cells in presence of both the cell traction mechanisms and the durotaxis mechanisms. During each time step, we used the Lemmon and Romer [43] model to calculate traction forces corresponding to current cell positions. Next, we started the finite element analysis from an undeformed matrix, calculating steady-state strains for the current traction forces. To simulate cell movement, which was biased by the local matrix strains using the durotaxis mechanism, we then applied one cell motility simulation time step, or Monte Carlo step (MCS; the MCS is the unit of time of our simulation; see Methods for detail and Discussion for an estimate of the real time corresponding to an MCS). After running the CPM for one MCS we again relaxed the matrix such that the next step started with an undeformed matrix. Thus we currently did not consider cell memory of substrate strains.
As Figure 3 and Video S1 demonstrate, in this model matrix stiffness affects both the morphology and motility of the simulated cells. On the most compliant substrate tested (0.5 kPa) the simulated cells contract and round up, whereas cells spread isotropically on the stiffest substrate tested (32 kPa). Overall, the cellular area increases with substrate stiffness (Figure 3 [39] and cardiac myocytes [36] in matrices of varying stiffness. The dependence of cell shapes on substrate stiffnesses is due to the transition from fluctuating to adherent pseudopods with increasing stiffness. Focal adhesions of cells on soft substrates all remain in the "fluctuating" state, irrespective of the local strains. On intermediate substrates, some pseudopods, due to increased traction, move to an extended state (mimicking a mature focal adhesion), generating more traction in this direction. Hence an initial stochastic elongation selfenhances in a feedback loop of increasing traction and strain stiffening. Such a self-enhancing cellelongation starting from an initial anisotropy in cell spreading has previously been suggested by Winer et al [44]. Extensions perpendicular to the long axis of an elongated cell do not occur since there is insufficient traction and the volume constraint is limiting. At matrices of high stiffness all pseudopods attempt to extend, mimicking the formation of static focal adhesion, until the volume constraint becomes limiting. This makes the cells spread more on stiff substrates than on soft substrates, with weaker volume constraints (lower values of λ) producing a stronger effect of substrate stiffness on cell shape and cell area ( Figure S1) .
We also measured the random motility of the cells by characterizing their dispersion coefficients, which we derived from the mean square displacements of the cells ( Figure S2; see section Morphometry for detail). The dispersion coefficients show biphasic behavior, with the highest motilities occurring at around 12 kPa (Figure 3 E). The biphasic dependence of the dispersion to substrate stiffness is in accordance with in vitro behavior of neutrophils [45], and smooth muscle cells [46]. Here it is typically thought to be due to a balance of adhesion and actin polymerization, or due to the interplay between focal adhesion dynamics and myosin-based contractility [45]. In our model, the effect is more likely due to the appearance of eccentric cell shapes at intermediate stiffnesses; as a result, only the tips of the cell generate sufficient strain in the matrix to extend pseudopods, producing more persistent motion than the round cells at stiff or soft substrates. It will be interesting to see if a similar relationship between cell shape and cell motility holds in vitro. Thus the model rules for cell traction and stretch guidance based on durotaxis and strain stiffening suffice to reproduce an experimentally plausible cellular response to matrix stiffness.

Mechanical cell-ECM feedback coordinates behavior of adjacent cells
Strains induced by endothelial cells on a compliant substrate with low concentrations of arginine-glycineaspartic acid(RGD)-containing nonapeptides can affect the behavior of adjacent cells [2]. On soft substrates (5.5 kPa or below) the cells reduced the motility of adjacent cells, whereas on stiff substrates (33 kPa) such an effect was not found. On substrates of intermediate stiffness (5.5 kPa), adjacent endothelial cells repeatedly attached and detached from one another, and cells moved more slowly in close vicinity of other cells, than when they were on their own. Because the extent to which cells could affect the motility of nearby cells depended on matrix compliancy, mechanical traction forces could act as a means for cell-cell communication [2]. To test if the simple strain-based mechanism represented in our model suffices for reproducing such mechanical cell-cell communication, we initiated the simulations with pairs of cells placed adjacent to one another at a distance of fourteen lattice sites corresponding to a distance of 35 µm, and ran a series of simulations on substrates of varying stiffness ( Figure 4 A and Video S2).
The cells behaved similar to the single cell simulations (Figure 3), with little cell-cell interactions at the lower and higher stiffness ranges. Consistent with previous observations [2], cell pairs on substrates of intermediate stiffness (12 kPa) dispersed more slowly than individual cells (paired two-sample t-test at 5000 MCS, p < 0.05 for 12 kPa), whereas individual cells and cell pairs dispersed at indistinguishable (p > 0.05) rates on stiff (14 kPa or more) or soft (10 kPa or below) substrates ( Figure 4, B-D) and Figure S3).
Also in agreement with the previous, experimental observations [2], on a simulated substrate of intermediate stiffness (12 kPa) the cells responded to the matrix strains induced by the adjacent cell by repeatedly touching each other, and separating again (Figure 4 E). The contact duration of cells on soft and stiff substrates, when they get close enough to each other, are typically longer than for intermediate substrates. This behavior is also similar to observations in vitro [2]. As one might expect that strongly adherent cells will not repeatedly touch and retract, but rather stay connected upon first contact, we investigated the effect of cell adhesion on these parameters ( Figure S4). Consistent with this intuition, for stronger adhesion, the contact count tends to be reduced and the contact durations tend to increase, but the overall trend holds: at intermediate matrix stiffnesses we continue to observe more frequent cell contacts than for more soft or more stiff matrices. Thus the observed pairwise cell behavior is primarily driven by durotaxis.
Mechanical strain can also coordinate the relative orientation of cells. Fibroblasts seeded on a compliant gel tend to align in a head-to-tail fashion along the orientation of mechanical strain [47]. Bischofs and Schwarz [48] proposed a computational model to explain this observation. Their model assumes that cells prefer the direction of maximal effective stiffness, where the cell has to do the least work to build up a force. This work is minimal between two aligned cells, because maximum strain stiffening occurs along the axis of contraction. Interestingly, visualization of our model results (Figure 1 C) suggested similar head-to-tail alignment of our model cells at around 12 kPa. To quantify cell alignment in our simulations, we measured the angle α between the lines l 1 and l 2 , defining the long axes of the cells and crossing the centers of mass of the cells (Figure 4 F). We classified the angles as acute (α < π/2; i.e. no alignment) or obtuse (α ≥ π/2; alignment). At matrix stiffnesses up to around 10 kPa, about one fourth of the angles α were obtuse, corresponding to the expected value for uncorrelated cell orientations. However, at 12 kPa and 14 kPa significantly more than a fourth of the angles α between the cell axes were obtuse (55/100 for 12 kPa, p < 1 × 10 −8 and 52/100 for 14 kPa, p < 1 × 10 −8 , binomial test), and for substrate compliancies of 8 to 16 kPa significantly more of the angles α were obtuse than for 4 kPa (p < 0.01 for 8 kPa, and p < 1 × 10 −12 for 10 kPa to 16 kPa; two-tailed Welch's t-test), suggesting that the mechanical coupling represented in our model causes cells to align in a head-to-tail fashion.

Mechanical cell-cell communication drives biologically-realistic collective cell behavior
After observing that the local, mechanical cell-ECM interactions assumed in our model sufficed for correctly reproducing many aspects of the behavior of individual endothelial cells on compliant matrices and of the mechanical communication of pairs of endothelial cells on compliant matrices, we asked what collective cell behavior the mechanical cell-cell coordination produced. When seeded subconfluently onto a compliant matrix (e.g., Matrigel), endothelial cells tend to organize into polygonal, vascular-like networks [5,6,49,50]. To mimic such endothelial cell cultures, we initialized our simulations with (approximately) 450 cells uniformly distributed over a lattice of 300×300 pixels (0.75×0.75 mm 2 ), corresponding to a cell density of 800 endothelial cells per mm 2 . In accordance with experimental observations on gels with low concentrations of collagen [6] or RGD-peptides [2], after 3000 MCS networks had not formed on soft matrices (0.5-4 kPa) or on stiff matrices (16-32 kPa) (  [5,6,49,50]. The optimal stiffness (≈ 10 kPa) for network formation is slightly lower than the stiffness of the substrate (≈ 12 kPa) on which single cells elongate the most (Figure 3 A). In comparison with a single cell, the collective pulling of a cell colony creates larger strains in the substrate. Consequently, the strain threshold inducing cell elongation is crossed at smaller substrate stiffness. that we compared the simulations of pairwise interactions with, in this experiment we used a 2.5 kPa gel functionalized with 5 µg/ml RGD peptide -a stiffness at which no network-formation is found in our simulations. Although we thus do not yet reach full quantitative agreement between model and experiment, note that network formation occurs at substrate stiffness of 10kPa on polyacrylamide matrices enriched with a low (1 µg/ml) concentration of collagen [6].
We next asked if the mechanical model could also reproduce sprouting from endothelial spheroids [7,8]. Video S5 and Figure 6 shows the results of simulations initiated with a two-dimensional spheroid of cells after 3000 MCS. On soft (0.5-8 kPa) and on stiff (32 kPa) matrices the spheroids stayed intact over the time course of the simulation. On matrices of intermediary stiffness (10-12 kPa) the spheroids formed distinct sprouts, visually resembling the formation of sprouts in in vitro endothelial spheroids [7,8]. On the 14 kPa and 16 kPa matrices the cells migrated away from the spheroid, with some cell alignment still visible for the 14 kPa matrices. Observation of a sprout protruding from a spheroid at 10 kPa suggests that a new sprout starts when one of the cells at the edge of the cluster protrudes and increases the strain in front of it. In a positive feedback loop via an increase in perceived stiffness the strain guides the protruding cell forward. The strain in its wake then guides the other cells along (Figure 6 C).

Discussion
In this paper we introduced a computational model of the in vitro collective behavior of endothelial cells seeded on compliant substrates. The model is based on the experimentally supported assumptions that (a) endothelial cells generate mechanical strains in the substrate [34,43], (b) they perceive a stiffening of the substate along the strain orientation, and (c) they extend preferentially on stiffer substrate [37]. Thus, in short, the assumptions are: cell traction, strain stiffening, and durotaxis. The model simulations showed that these assumptions suffice to reproduce, in silico, experimentally observed behavior of endothelial cells at three higher level spatial scales: the single cell level, cell pairs, and the collective behavior of endothelial cells. In accordance with experimental observation [36,39], the simulated cells spread out on stiff matrices, they contracted on soft matrices, and elongated on matrices of intermediate stiffness (  [2]. Also, in agreement with experimental observations of fibroblasts on compliant substrates [47] and previous model studies [48] the cells repositioned into an aligned, head-to-tail orientation (Figure 4 F). The model simulations further suggest that these pairwise cell-cell interactions suffice for vascular-like network formation in vitro ( Figure 5) and sprouting of endothelial spheroids ( Figure 6).
The correlation between pairwise cell-cell interactions and collective cell behavior observed in our computational model parallels observations in vitro. Cells elongate due to positive feedback between stretch-guided extension and cell traction, as previously suggested by Winer et al. [44]. Elongated and spindle-shaped cells are considered indicative of future cell network assembly [6]. Our model suggests that the elongated cell shapes produce oriented strains in the matrix, via which cells sense one another at a distance. In this way new connections are continuously formed over "strain bridges" (see, e.g., Figure 5 C,D and Video S4), while other cellular connections break, producing dynamically stable networks as illustrated in Video S3. Such dynamic network restructuring was also observed during early embryonic development of the quail embryo [51] and in bovine aortic endothelial cell cultures ( Figure 5 D and [6]), but not in human umbilical vein endothelial cell cultures [23,50]. Also in agreement with experimental results, the collective behavior predicted by our model strongly depends on substrate stiffness. The strongest interaction between cell pairs is found on substrates of intermediate stiffness, enabling network formation [2], whereas network assembly does not occur on stiffer or on softer substrates [6].
These agreements with experimental results are encouraging, but our model also lacks a number of properties of in vitro angiogenesis that pinpoint key components still missing from our description. We compared the simulation of pairwise cell-cell interactions with previous experiments conducted on polyacrylamide gels, functionalized with RGD ligands [2], which have linear elastic behavior for small deformations [52][53][54]. Strain-stiffening of polyacrylamide gels has been reported for deformations over 2 µm [55]. Thus with pixels in our model measuring 2.5 µm × 2.5 µm, strain-stiffening seems a reasonable assumption. Nevertheless, a possible alternative interpretation of the cell pair simulations is that the increased tension generated in pseudopods pulling on the matrix leads to a higher probability of focal adhesion maturation [41,42]. A further issue is that in our simulations, single cells dispersed somewhat more quickly on soft gels than on stiff gels (Figure 3 E and Figure S2). This model behavior contradicts experimental observations that endothelial cells move fastest on stiff substrates [2]. Another open issue concerns the time scales of our simulations. In the present paper time we use the Monte Carlo step as a (computational) unit of time. To estimate the actual time corresponding to 1 MCS, we scale the single cell dispersion coefficients shown in Figure 3 E to experimental dispersion coefficients of bovine endothelial cells on compliant substrates in vitro [2]. Reported dispersion coefficients of endothelial cells range from around 1 µm 2 /min (on substrates of 500 Pa) to around 10 µm 2 /min (on substrates of 5500 Pa) (as derived from the MSDs in Figure 3a,c in [2] and based on MSD(t) = 4Dt; cf. Eq. 13). The dispersion coefficients of single cells in our simulations are in the range of 0.03 − 0.08 µm 2 /MCS (Figure 3), assuming pixels of 2.5 × 2.5 µm 2 . Thus, based on fitting of single cell dispersion rates, the estimated length of 1 MCS is 0.5 to 3 seconds. The typical time scale of a vascular network formation simulation is around 3000 MCS ( Figure 5), i.e., 12.5 min to 2.5 hr for these time scale estimates. In experiments, network formation takes longer, around 24 hr. Thus in our current model the time scales of cell dispersion and network formation do not match exactly. A possible reason of this discrepancy is the short persistent length of cell motility in standard cellular Potts models. To better match the time scales of single cells and collective cell behavior in our model, in our future work we will increase the persistence length of the endothelial cells by using the available cellular Potts methodology [56][57][58], or model the subcellular mechanisms of cell motility in more detail, e.g. by including mean-field models of actin polymerization [59,60]. A further open issue is the interaction between substrate mechanics and cell-substrate adhesivity. Although the model correctly predicts the absence of network formation on stiff substrates, it cannot yet explain the observation that reducing the substrate adhesivity of the endothelial cells rescues network formation on stiff substrates [6]. On compliant gels endothelial cells must secrete fibronectin to form stable networks, whereas fibronectin polymerization inhibitors elicit spindle-like cellular phenotypes associated with network formation on stiff matrices, under conditions where networks do not normally form [6]. To explain these observations, straightforward future extensions of the model will include a more detailed description of cell-substrate adhesion, combined with models of ECM secretion and proteolysis [13,25,27,61].
The current model also assumes a uniform density (i.e., the infinitesimal strain assumption) and thickness of the extracellular matrix, whereas under some culture conditions the endothelial cells have been reported to pull the extracellular matrix underneath them [62], producing gradient in matrix density and/or thickness. To describe the role of viscous deformations of the extracellular matrix in morphogenesis, Oster and Murray [17,18] developed a continuum mechanical model of pattern formation in mesenchymal tissues. Their model assumed (a) that cells exert contractile forces onto the surrounding extracellular matrix, that will (b) locally deform the ECM, resulting in passive displacements of cells along with the ECM, and (c) produce density gradients in the ECM along which cells move actively due to haptotaxis. These mechanisms together produce periodic cell density patterns. Manoussaki et al. [15] and Namy et al. [19] applied this work to investigate mechanical cell-ECM interactions during angiogenesis, and demonstrated that the mechanism can produce vascular-like network patterns. In their model they also included an anisotropic diffusion term to simulate preferential movement along the local strain-direction, but the term was neither necessary nor sufficient for network formation. This finding contradicts our model in which strain-induced sprouting is the driving force of network formation and sprouting. The two models represent the two extremes of network formation on visco-elastic matrices. Here, the Manoussaki et al. [15] and Namy et al. [19] models represent patterning on viscous matrices, in which cellular traction forces pull the matrix together while inducing little strain or stress. Our model would represent elastic materials, in which pulling forces induce local strains. Future extensions of the model will include matrix remodelling (e.g., by assuming a matrix thickness field) allowing us to study the full range of viscoelastic matrices.
Apart from these biological issues, we made several mathematical simplifications that we will improve upon in future models of cell-ECM interactions. In the current model, for mathematical simplicity, we assumed that after each Monte Carlo step the matrix was undeformed again. Thus we currently did not consider cell memory of substrate strains. Further developments of the model presented here will improve on this issue, because actin filament dynamics are typically influenced by the past evolution of substrate deformations, e.g., due to reorientation of matrix fibers [62]. For computational efficiency, we assumed linearly elastic materials and infinitesimal strain in the finite element simulations, and mimicked durotaxis via a perceived strain-stiffening (Eq. 9) where cells perceive increased ECM stiffness due to local strain. In our ongoing work we are interfacing the open source package FEBio (http://febio.org) with the cellular Potts package CompuCell3D (http://compucell3D.org). This will allow us to run our model with any ECM material available to users of FEBio, including strain-stiffening materials. Using an actual strain stiffening material may lead to longer-range interactions between cells, because locally stiffer regions may channel the stress between the cells [63]. A further technical limitation of our model is that we currently only run two-dimensional simulations, representing cells moving on top of a two-dimensional culture system. The ongoing interfacing of FEBio and CompuCell3D will pave the way for modeling cell-ECM interactions in three-dimensional tissue cultures. We also plan to model fibrous extracellular matrix materials in more detail.
A quite puzzling aspect of vascular network formation and spheroid sprouting is that so many alternative, often equally plausible computational models can explain it (reviewed in [12]). Including the present model, there are at least three alternative computational models based on mechanical cell-ECM interactions [15,16,19,31,64], a series of models assuming chemoattraction between endothelial cells [20,21,23,24,65,66] and extensions thereof [25,27,67], and models explaining network formation in absence of chemical or mechanical fields [29,30,32]. Each of the models explains one aspect of vascular network formation or a response to an experimental treatment that the other models cannot explain, e.g. the relation between spindle-shaped cell phenotypes and network formation [23,30], the requirement of VE-cadherin signaling for network formation and sprouting [24,29], the binding and release of growth factors from the ECM [25,26], the role of mechanical ECM restructuring and haptotaxis [15,19,31], the response of vascular networks to toxins [27], or the role of intracellular Ca 2+ signaling [57]. Among these alternative models, we must now experimentally falsify incorrect mechanisms, and fine-tune and possibly combine the remaining models to arrive at a more complete understanding of the mechanisms of angiogenesis. To this end, we are currently quantitatively comparing the kinetics of patterns produced by chemotaxis-based, traction-based, and cell-elongation based models with the kinetics of in vitro networks [23,50]. The resulting, more complete model would likely contain aspects of each of the available computational models and assist in explaining the conflicting results obtained from the available experimental systems, culture conditions, and in silico models of angiogenesis.

Methods
To model the biomechanical interactions between endothelial cells and compliant matrices, we developed a hybrid of the cellular Potts model (CPM) [68,69] to represent the stochastic motility of the endothelial cells, and a mechanical model based on the finite element method (FEM) [70] of the compliant extracellular matrix. Related CPM-FEM models were proposed for the simulation of load-induced bone remodeling [71,72], and recently a related approach was proposed in a model study of cell alignment [73]. A documented simulation code is provided as part of the Supporting Information (Supporting Text S1 and Code S1) and a detailed list of parameter values is given in Table S1.

Cellular Potts model
The CPM represents cells on a regular square lattice, with one biological cell covering a cluster of connected lattice sites. To mimic random cell motility, the CPM iteratively expands and retracts the boundaries of the cells, depending on the passive forces acting on them and on the active forces exerted by the cells themselves. These are summarized in a balance of forces, represented by the Hamiltonian, The first term is an (approximate) volume constraint, with a(σ) the actual volume of the cells, A(σ), a resting volume, and λ an elasticity parameter that regulates the permitted fluctuation around the resting volume. In contrast with the original formulation of the CPM [68], the deviation of the cell from its target volume is taken relative to the target volume, by analogy with the (non-dimensional) engineering strain. Alternative, similar volume constraints can be chosen [67]. We use a value A(σ) = 50 for all cells; the medium does not have a volume constraint. The second term represents cell-cell and cellmedium adhesion, where J(σ( x), σ( x )) ≥ 0 is the contact cost between two neighboring pixels, and δ, the Kronecker delta. Throughout the manuscript we use neutral cell-cell adhesion settings; J(σ( x), σ( x )) = 2.5 at cell-cell interfaces, and J(σ( x), 0) = 1.25 at cell-medium interfaces, with σ( x) > 0 and σ( x ) > 0. In other words, cells have no preference for adhering to other cells or the medium. For these neutral cell adhesion parameter settings, cells will still adhere weakly to one another (a remedy for this effect was proposed in [74]). Additional terms in the Hamiltonian represent the cells' responses to ECM mechanics, and will be described in more detail below. The CPM iteratively selects a random lattice site x and attempts to copy its state, σ( x ), into a randomly selected adjacent lattice site x. To reflect the physical, "passive" behavioral response of the cells to their environment, the copy step is always accepted if it decreases the Hamiltonian. To account for the active random motility of biological cells, we allow for energetically unfavorable cell moves, by accepting copies that increase the Hamiltonian with Boltzmann probability, where ∆H is the change in H if the copying were to occur, and T > 0 parameterizes the intrinsic cell motility. It represents the extent to which the active cell motility can overcome the reactive forces (e.g. volume constraint or adhesions) in the environment. We assume that all cells keep the same motility and thus set T to be constant throughout the simulations. During one Monte Carlo step (MCS), we perform n copy attempts, with n equal to the number of sites in the lattice. To prevent cells from splitting up into two or more disconnected patches, we use a connectivity constraint that rejects a spin flip σ( x ) → x if it would break apart the retracting cell σ( x).

Model of Compliant Substrate based on Finite Element Method
A two-dimensional model describes the compliant substrate on which the cells move. Deformations are calculated using the finite element method (FEM; reviewed in [70]). The FEM represents the substrate as a lattice of finite elements, e, with each element corresponding to a pixel of the CPM. To obtain the finite element equations, the weak formulation (associated with the total potential energy) of the governing equations of the displacement u of the substrate is set up, in order to obtain the finite element equations, with stiffness matrix K, displacement u, and forces f . The vector u = [u x1 , u y1 , u x2 , ...u xn , u yn ] T contains the displacements of all nodes, which are the unknowns that the FEM calculates based on the active forces exerted onto the material, presented in f . In this paper f only consists of traction forces that the cells apply onto the ECM, unless stated otherwise. In a two-dimensional analysis the forces f are divided by the thickness they are working on. For this we assume an effective substrate thickness t = 10 µm. We impose boundary conditions of u = 0 at the boundary of the CPM grid, this means that the substrate is fixed along the boundaries.
To a first approximation, in this work we consider an isotropic, uniform, linearly elastic substrate [48,75] and we apply infinitesimal strain theory: We assume that material properties, including local density and stiffness are unchanged by deformations. The global stiffness matrix K is assembled from the element stiffness matrices K e (see Supporting Text S1 and [70]), which describe the relation between nodes of each element, e, where B-the conventional strain-displacement matrix for a four-noded quadrilateral element (see Supporting Text S1 and [70])-relates the node displacements u e to the local strains, as, The strain vector is a column notation of the strain tensor and D is the material property matrix. Assuming plane stress conditions, where E is the material's Young's modulus, and ν is Poisson's ratio. Throughout this study, we use a Poisson's ratio ν = 0.45 and Young's moduli ranging from E = 0.5 kPa to E = 32 kPa, which are plausible values for most cell culture substrates [48,53,76]. For more details of the derivation of Eq. 3, and the entries in B, see Supporting Text S1 and [70]. As a reference configuration for the displacements we used an unstretched substrate, u = 0. Thus, after each Monte Carlo step (during which the cells positions and shapes have changed) the substrate is assumed to be undeformed, such that the stiffness matrix, K, is constant in time. This prevents expensive calculations that would be necessary for recalculating K in each iteration. Although the previous displacements do not influence the new deformation of the substrate, they are used as an initial guess for solving Ku = f , in order to reduce the number of iterations necessary to converge to the FEM solution.

Mechanical cell-substrate coupling
To simulate cell-substrate feedback we alternate the cellular Potts model (CPM) steps with a simulation of the substrate deformations using the finite element method. We assume that cells apply a cell-shape dependent traction on the ECM and the cells respond to the resulting ECM strains by adjusting their cell shape. Using the CPM grid as the finite element mesh, the pixels of the CPM become four-node square elements in the FE-mesh. Adopting the model by Lemmon & Romer [43], we assume that each node i covered by a CPM cell pulls on all other nodes j in the same cell, at a force proportional to distance d i,j . The resultant force F i on node i then becomes, where ∆x is the lattice spacing and µ gives the tension per unit length. This parameter has been scaled to µ = 0.01 nN/µm, such that the total cell traction corresponds to experimentally reported values [77]. The resultant forces point towards the cell centroid, and are proportional to the distance from it ( Figure 2). In this way a CPM configuration yields a traction force F , which are collected in the forces f for the finite element calculation. To calculate the resulting ECM strains, we solve Ku = f for the node displacements u with a preconditioned conjugate gradient (PCG) solver [78], and derive the local strains using Eq. 5. The reference configuration for the displacements is an unstretched substrate, u = 0. After a sufficiently accurate solution for the FEM equations has been obtained by the PCG solver, we run a Monte Carlo step of the CPM. After each MCS, which changes cell positions, the substrate is assumed to be undeformed again, for the sake of simplicity. Thus, the stiffness matrix, K, is constant in time.
We assume durotaxis, i.e., the CPM cells preferentially extend pseudopods on matrices of higher stiffness (e.g., because of strain stiffening). By analogy with chemotaxis algorithms [79] at the time of copying we add the following durotaxis term to ∆H in response to the strain-and orientation-dependent ECM stiffness E, with g( x, x ) = 1 for extensions and g( x, x ) = −1 for retractions, λ durotaxis is a parameter, v m = x − x , a unit vector giving the copy direction, and 1 and 2 , and v 1 and v 2 eigenvalues and eigenvectors of representing the principal strains and strain orientation. We use the strain ( x) in the target pixel when considering an extension, and for retractions we use the strain in the source pixel, ( x ). Thus we consider the strain in the ECM adjacent to the pseudopod. The sigmoid h(E) = 1/(1+exp(−β(E −E θ ))), with threshold stiffness E θ , and β, the steepness of the sigmoid, mimics maturation of focal adhesions under the influence of tension [41]. The tension in focal adhesions will increase with higher local matrix stiffness, E, because the matrix will deform less easily. The sigmoid function starts at zero, goes up when there is sufficient stiffness, and eventually reaches a maximum. This means that a certain level of stiffness is needed to cause a cell to spread. Alternative forms of h(E) can be used: For an overview see Figure S5. Due to limitations of our current finite element code and for reasons of computational efficiency, we assumed a linearly elastic, isotropic material in the FEM, thus precluding explicit strain stiffening effects in the FEM calculations. Instead, we implemented the effect of strain-stiffening in the cell response, where cells perceive increased ECM stiffness as a function of the principal strains 1 and 2 , where E 0 sets a base stiffness for the substrate, and st is a stiffening parameter. The indicator function 1 >0 = {1, > 0; 0, ≤ 0} indicates that strain stiffening of the substrate only occurs for substrate extensions ( ≥ 0); compression ( < 0) does not stiffen or soften the substrate.

Morphometry
To characterize the random motility of single cells and cell pairs, we measured the cells' mean square displacement, with C(S, t), the centroid of cell S at Monte Carlo step ("time") t, given by with C(S, t), the set of coordinates of the lattice sites comprising cell S at MCS t, and x = {x 1 , x 2 }. The MSD is a reliable measure of random motility [80] and it can be directly compared with experimental data (e.g., [2]).
The dispersion coefficient, defined as is derived from the slope of the MSD, and is used as a measure of the motility of random walkers. The length, orientation and eccentricity of cells were estimated from the inertia tensors I(S) of the cells, defined as [81], Assuming

Endothelial Cell Culture
Bovine aortic endothelial cells (BAECs) (VEC Technologies, Rensselaer, NY) were cultured through passage 12. Cells were kept at 37 • C and 5% CO 2 and fed every other day with Medium 199 (Invitrogen, Carlsbad, CA) supplemented with 10% Fetal Clone III (HyClone, Logan, UT), 1% MEM amino acids (Invitrogen), 1% MEM vitamins (Medtech, Manassas, VA), and 1% penicillin-streptomyocin (Invitrogen). Polyacrylamide hydrogels were synthesized as previously described [6]. Briefly, a gel mixture was prepared from MilliQ water, HEPES, TEMED (Bio-Rad, Hercules, CA) and a 5%:0.1% ratio of acrylamide to bisacrylamide (Bio-Rad) to generate substrates with a Young's modulus of 2,500 Pascals. Polymerization was initiated by the addition of N-6-((acryloyl)amido)hexanoic acid (synthesized according to Pless et al. [82]) and ammonium persulfate (Bio-Rad) into the gel mixture. Following polymerization, gels were incubated with 5 µg/ml RGD peptide (GCGYGRGDSPG) (Genscript), followed by ethanolamine (Sigma). Gels were stored in PBS overnight. Hydrogels were sterilized with ultraviolet light before cell culture. A T-75 flask with a confluent BAEC monolayer was seeded onto the hydrogels at 350,000 cells per gel (approximately 1,375 cells per mm 2 ). The gels were maintained at 37 • C and 5% CO 2 for three days prior to imaging. After replenishing with fresh complete media, the cells on hydrogels were visualized with a Zeiss Axio Observer.Z1 inverted spinning disc microscope with a Hamamatsu ORCA-R 2 digital camera. Images were captured every 30 minutes for 24 hours.   Supporting Information Table S1. Parameter settings of the simulation model.     Video S1. Behavior in silico of a single cell on substrates of 4 kPa, 12 kPa, and 32 kPa, for a duration of 500 MCS per simulation. Parameter settings as in Figure 3.
Video S2. Pairwise cell-cell interactions in silico on substrates of 4 kPa, 12 kPa, and 32 kPa, for a duration of 500 MCS per simulation. Parameter settings as in Figure 4.
Video S3. Network formation in silico on a substrate of 10kPa, for a duration of 3000 MCS. Video represents a 0.75 × 0.75 mm 2 area (300 × 300 pixels) initiated with 450 cells. Parameter settings are as in Figure 5.
Video S4. Network formation of bovine aortic endothelial cells on a 2.5 kPa polyacrylamide gel functionalized with RGD-peptide. Time lapse images were captured in 30 minute intervals over an 8 hour time period. Image size as in Figure 5 D.
Video S5. Sprouting in silico from a spheroid on a substrate of 10kPa, for a duration of 3000 MCS. Video represents a 0.75 × 0.75 mm 2 area (300 × 300 pixels) initiated with 450 cells. Parameter settings are as in Figure 6.
Supporting Text S1. Documentation of C and Matlab code used for the simulations, including a detailed description of the finite-element model.
Protocol S1. C and Matlab source code used for the simulations.