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.

Introduction

How the behavior of cells in a multicellular organism is coordinated to form structured tissues, organs and whole organisms, is a central question in developmental biology. Keys to answering this question are chemical and mechanical cell-cell communication and the biophysics of self-organization. Cells exchange information by means of diffusing molecular signals, and by membrane-bound molecular signals for which direct cell-cell contact is required. In general, these developmental signals are short-lived and move over short distances. The extracellular matrix (ECM), the jelly or hard materials that cells secrete, provides the micro-environment the cells live in. Apart from its supportive function, the ECM mediates molecular [1] and biomechanical [2] signals between cells. Mechanical signals, in the form of tissue strains and stresses to which cells respond [3], can act over long distances and integrate mechanical information over the whole tissue [4], and also mediate short-range, mechanical cell-cell communication [2]. How such mechanical cell-cell communication via the ECM can coordinate the self-organization of cells into tissues is still poorly understood. Here we propose a cell-based model of endothelial cell motility on compliant matrices to address this problem.

A widely used approach to study the role of cell-ECM interactions in coordinating collective cell behavior is to isolate cells (e.g., endothelial cells isolate from bovine aortae or from human umbilical cords or foreskins) and culture them on top of or inside an artificial or natural ECM (e.g., Matrigel). This makes it possible to study the intrinsic ability of cells to form tissues in absence of potential organizing signals or pre-patterns from adjacent tissues. A problem particularly well-studied in cell cultures is the ability of endothelial cells to form blood vessel-like structures, including the formation of vascular-like networks from dispersed cells and the sprouting of spheroids. To this end, cell cultures can be initialized with a dispersion of endothelial cells on top of an ECM material (e.g., Matrigel, collagen, or fibrin) [5], [6], with endothelial spheroids embedded within the ECM [7], [8], or with confluent endothelial monolayers [9]–[11]. Although the conditions required for vascular-like development in these in vitro culture systems are well established, the mechanisms driving pattern formation of endothelial cells are heavily debated, and a wide range of plausible mechanisms has been proposed in the form of mathematical and computational models reproducing aspects of angiogenesis (reviewed in [12]–[14]).

Typical ingredients of network formation models are (a) an attractive force between endothelial cells, which is (b) proportional to the cell density, and (c) inhibited or attenuated at higher cellular densities. The attractive force can be due to mechanical traction or due to chemotaxis. Manoussaki, Murray, and coworkers [15], [16] proposed a mechanical model of angiogenic network formation, based on the Oster and Murray [17], [18] continuum mechanics theory of morphogenesis. In their model, endothelial cells exert a uniform traction force on the ECM, dragging the ECM and the associated endothelial cells towards them. The traction forces saturated at a maximum cell density. Namy and coworkers[19] replaced the endothelial cells' passive motion along with the ECM for active cell motility via haptotaxis, in which cells move actively towards higher concentrations of the ECM. Both models also included a strain-biased random walk term for the endothelial cells, but they found that it had little effect on network formation; the mechanism was dominated by cell aggregation. In their model based on chemotaxis, Preziosi and coworkers [20], [21] assumed that cells attract one another via the secreted chemoattractant VEGF. Due to diffusion and first-order degradation, the chemoattractant forms exponential gradients around cells leading to cell aggregation in much the same way as that assumed in the Manoussaki and Namy models. These chemotaxis-based hypotheses formed the basis for a series of cell-based models based on the cellular Potts model (CPM). Assuming chemotactic cell-cell attraction, and a biologically-plausible overdamped cell motility, the cells in these CPM models form round aggregates, in accordance with the Keller-Segel model of cell aggregation [22]. Additional assumptions, including an elongated cell shape [23] or contact 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]–[27] and a range of additional secreted and exogenous growth factors [27], and studied the ability of the contact-inhibition 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.

Results

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 (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 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].

Figure 1

Simulated cellular responses to static strains.

Cells do not generate traction forces in this figure. (A) Cell length as a function of the durotaxis parameter, , on a substrate stretched along the vertical axis. (B) Cell orientation as a function of the stretch orientation (simulated with ). Error bars show standard deviation for . Insets show five simulations per value tested.

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 ,. 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].

Figure 2

Visualization of simulated traction forces (black arrows) and resulting matrix strains (blue line segments) generated in the proposed hybrid cellular Potts and finite element simulation model.

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 B). On matrices of intermediate stiffnesses (around 12 kPa) the cells elongate, as reflected by measurements of the cell length (Figure 3 C) and eccentricity (Figure 3 D) that both have maximum values at around 12 kPa. Such a biphasic dependence of cellular morphology on the stiffness of the ECM mimics the behavior of endothelial cells [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 self-enhances in a feedback loop of increasing traction and strain stiffening. Such a self-enhancing cell-elongation 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).

Figure 3

Simulated individual cell responses to mechanical cell-ECM feedback.

(A) Single cells on substrates of varying stiffness after 100 MCS. Line pieces indicate strain magnitude and orientation. (B) cell area () of cells; (C) cell length (length of major axis if the cell is seen as an ellipse) as a function of substrate stiffness (D) cell eccentricity (, with and the lengths of the cell's major and minor semi-axes) as a function of stiffness. Mean and standard deviation shown for in panels B-D. (E) Dispersion coefficients of individual, simulated cells, derived from a linear fit on the mean square displacements (Figure S2); . Error bars indicate 95% confidence intervals of linear fits.

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-glycine-aspartic 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).

Figure 4

Simulated cell-cell interactions on substrates of varying stiffnesses.

(A) Visualization of cell shapes and substrate strains in absence of external strain. Line pieces indicate strain magnitude and orientation. (B-D) Mean square displacement of individual cells (blue errorbars) and cell pairs (red errorbars) on simulated substrates. (B) 4 kPa; (C) 12 kPa; (D) 32 kPa. Error bars in panels B to D indicate standard deviation for . (E) Number of cell-cell contacts made over 500 MCS between two simulated cells initiated at a distance of fourteen lattice sites from each other. Error bars show standard deviation over simulations (F) Quantification of head-to-tail alignment of cells. An obtuse angle between the two cells' long axes indicates that cells are oriented head-to-tail. Plotted is the fraction of Monte Carlo steps over MCS 20-500 that the two cells are aligned head-to-tail. Shown are means and standard deviations over 100 independent simulations on a field of 0.25×0.25 (100×100 pixels). Insets: examples of acute (left) and obtuse (right) cell configurations.

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 and , defining the long axes of the cells (Figure 4 F). We classified the angles as acute (; i.e. no alignment) or obtuse (; 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.75mm2), corresponding to a cell density of 800 endothelial cells per mm2. 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) (Figure 5 A): The cells tended to form small clusters (Figure 5 A). Interestingly, on matrices of intermediate stiffness after around 300 MCS the cells organized into chains (8 kPa) or network-like structures (10 kPa and 12 kPa) similar to vascular network-like structures observed in endothelial cell cultures [5], [6], [49], [50]. The optimal stiffness (≈10kPa) for network formation is slightly lower than the stiffness of the substrate (≈12kPa) 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.

Figure 5

Simulated network formation assay.

(A) Simulated collective cell behavior on substrates of varying stiffness, with a uniformly distributed initiated configuration of cells. (B) Time lapse showing the development of a polygonal network on a 10kPa substrate (time in MCS). Panels A and B represent a 0.75×0.75 area ( pixels) initiated with 450 cells. (C) Close-up of simulated network formation on a 10 kPa substrate, showing the reconnection of two sprouts. Time in MCS. (D) Time lapse imaging of bovine aortic endothelial cells seeded onto a 2.5 kPa polyacrylamide gel functionalized with RGD-peptide. Arrows indicate cells that join together and elongate into a network. Time scale is in hours. Scale bar is 50 µm.

Figure 5 B and Video S3 show a time-lapse of the development of a network configuration on a substrate of 10kPa. The cells organized into a network structure within a few hundred MCS. The network was dynamically stable, with minor remodeling events taking place, including closure and bridging of lacunae. Figure 5 C shows such a bridging event in detail. In an existing lacuna (1800 MCS) stretch lines bridged the lacuna, and connected two groups of cells penetrating the lacuna (1980 MCS). The cells preferentially followed the path formed by these stretch lines (2150 MCS) and reached the other side of the lacuna by 2400 MCS. Such bridging events visually resemble sprouting in bovine endothelial cell cultures on compliant matrices (Figure 5 D, Video S4, and [6]). To stay close to the experimental conditions used for the observations of pairwise endothelial cell-cell interaction on compliant substrates [2] 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).

Simulated spheroid assay.

(A) Collective behavior in a simulation initiated with a two-dimensional "spheroid" of cells, on substrates of varying stiffness. (B) Time lapse showing a sprouting spheroid on a 10kPa substrate. Time in MCS. Panels A and B represent a 0.75×0.75 mm2 area (300×300 pixels) initiated with a spheroid consisting of 113 cells; (C) Close-up of sprouting on a 10 kPa substrate. Time in MCS. Black line pieces indicate strain magnitude and orientation.

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 (Figure 3). The same assumptions also sufficed to reproduce experimentally observed pairwise cell-cell coordination. On matrices of intermediate stiffness, endothelial cells slowed down each other (Figure 4 B) and repeatedly touched and retracted from each other (Figure 4 E and Video S2), in agreement with in vitro observations of bovine aortic endothelial cells on acrylamide gels [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]–[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 (on substrates of ) to around (on substrates of ) (as derived from the MSDs in Figure 3a,c in [2] and based on ; cf. Eq. 13). The dispersion coefficients of single cells in our simulations are in the range of (Figure 3), assuming pixels of . 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., to 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]–[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 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 the actual volume of the cells, , 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 for all cells; the medium does not have a volume constraint. The second term represents cell-cell and cell-medium adhesion, where is the contact cost between two neighboring pixels, and , the Kronecker delta. Throughout the manuscript we use neutral cell-cell adhesion settings; at cell-cell interfaces, and at cell-medium interfaces, with and . 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 and attempts to copy its state, , into a randomly selected adjacent lattice site . 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 is the change in H if the copying were to occur, and 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 to be constant throughout the simulations. During one Monte Carlo step (MCS), we perform copy attempts, with 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 if it would break apart the retracting cell .

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, , 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 of the substrate is set up, in order to obtain the finite element equations,

with stiffness matrix , displacement , and forces . The vector 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 . In this paper only consists of traction forces that the cells apply onto the ECM, unless stated otherwise. In a two-dimensional analysis the forces are divided by the thickness they are working on. For this we assume an effective substrate thickness . We impose boundary conditions of 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 is assembled from the element stiffness matrices (see Supporting Text S1 and [70]), which describe the relation between nodes of each element, ,

where —the conventional strain-displacement matrix for a four-noded quadrilateral element (see Supporting Text S1 and [70])—relates the node displacements to the local strains, as,

The strain vector is a column notation of the strain tensor and is the material property matrix. Assuming plane stress conditions,

where is the material's Young's modulus, and is Poisson's ratio. Throughout this study, we use a Poisson's ratio and Young's moduli ranging from to , 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 , see Supporting Text S1 and [70].

As a reference configuration for the displacements we used an unstretched substrate, . 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, , is constant in time. This prevents expensive calculations that would be necessary for recalculating 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 , 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 covered by a CPM cell pulls on all other nodes in the same cell, at a force proportional to distance . The resultant force on node then becomes,

where is the lattice spacing and gives the tension per unit length. This parameter has been scaled to , 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 , which are collected in the forces for the finite element calculation. To calculate the resulting ECM strains, we solve for the node displacements 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, . 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, , 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 in response to the strain- and orientation-dependent ECM stiffness ,

with for extensions and for retractions, is a parameter, , a unit vector giving the copy direction, and and , and and eigenvalues and eigenvectors of representing the principal strains and strain orientation. We use the strain in the target pixel when considering an extension, and for retractions we use the strain in the source pixel, . Thus we consider the strain in the ECM adjacent to the pseudopod. The sigmoid , with threshold stiffness , 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, , 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 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 and ,

where sets a base stiffness for the substrate, and is a stiffening parameter. The indicator function indicates that strain stiffening of the substrate only occurs for substrate extensions (); compression () 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 , the centroid of cell at Monte Carlo step (“time”) , given by

with , the set of coordinates of the lattice sites comprising cell at MCS ,

and . 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 of the cells, defined as [81],

Assuming cells are approximately ellipse-shaped, the length of cell is approximated as , with the largest eigenvalue of . The eccentricity of a cell is defined using the eigenvalues of the inertia tensor as , where are the eigenvalues of . An eccentricity close to zero corresponds to roughly circular cells and cells with an eccentricity close to unity are more elongated. The orientation of the cell is given by the eigenvectors of the inertia tensor .

Endothelial cell culture

Bovine aortic endothelial cells (BAECs) (VEC Technologies, Rensselaer, NY) were cultured through passage 12. Cells were kept at and 5% 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 bis-acrylamide (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 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 mm2). The gels were maintained at and 5% 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-R2 digital camera. Images were captured every 30 minutes for 24 hours.

Simulated responses of individual cells to mechanical cell-ECM feedback as a function of the values of the volume restriction, . Columns: area (left), cell length (middle) and eccentricity (right). Mean and standard deviation shown for after 500 MCS on simulated substrates of stiffness varying from 0.5 kPa to 32 kPa.

(PDF)

Click here for additional data file.

Mean square displacements of individual cells on simulated substrates of stiffness varying from 0.5 kPa to 32 kPa. Mean square displacement shown over cells.

(PDF)

Click here for additional data file.

Mean square displacement of individual cells (blue errorbars) and cell pairs (red errorbars) on simulated substrates of stiffness varying from 0.5 kPa to 32 kPa. Error bars indicate standard deviation for .

(PDF)

Click here for additional data file.

Number of cell-cell contacts made over 500 MCS (left column) and contact duration (right column) over 500 MCS between two simulated cells initiated at a distance of fourteen lattice sites from each other on simulated substrates of stiffness varying from 0.5 kPa to 32 kPa, for intercellular contact energies varying from

(PDF)

Click here for additional data file.

Effect of form of model function Standard, sigmoid function, as used in main text, with , , and . (B) Saturated function, , with and . (C) Piecewise linear function, , with , , and . (D) Gaussian function, , with and , . Insets show typical cell shape for regions indicated by red bars.

(PDF)

Click here for additional data file.

C and Matlab source code used for the simulations.

(ZIP)

Click here for additional data file.

Parameter settings of the simulation model.

(PDF)

Click here for additional data file.

Documentation of C and Matlab code used for the simulations, including a detailed description of the finite-element model.

(PDF)

Click here for additional data file.

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.

(MOV)

Click here for additional data file.

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.

(MOV)

Click here for additional data file.

Network formation in silico on a substrate of 10kPa, for a duration of 3000 MCS. Video represents a area ( pixels) initiated with 450 cells. Parameter settings are as in Figure 5.

(MOV)

Click here for additional data file.

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.

(MOV)

Click here for additional data file.

Sprouting in silico from a spheroid on a substrate of 10kPa, for a duration of 3000 MCS. Video represents a area ( pixels) initiated with 450 cells. Parameter settings are as in Figure 6.

(MOV)

Click here for additional data file.