Optimal mechanical interactions direct multicellular network formation on elastic substrates

To model cell network formation, we incorporate substrate-mediated cell mechanical interactions into an agent-based model for cell motility (45). This captures the dynamic rearrangements of cells into favorable configurations. In our agent-based approach (46, 47), summarized in Fig. 1B, we consider a system of $N$ particles, each a disk of diameter $d$. Depending on the context, each disk could model a cell or its constituent parts, and their motion represents both cell migration as well as cell spreading or shape change dynamics. Details of the cell shape are not included in this minimal model. These model cells self-organize according to substrate friction–dominated overdamped dynamics that depend on intercell interactions as well as individual cell stochastic movements described by an effective diffusion. The model incorporates both short-range, steric and long-range, substrate-mediated elastic interactions between cells and is detailed in the Materials and Methods.

The ubiquitous traction force pattern generated by a single polarized cell with a long axis $\mathbf{a}$ and exerting a typical force $\mathbf{F}$ at its adhesions can be modeled as a force dipole, $P_{\mathit{ij}}=F_i{a}_j$ (Fig. 1C). Note that the cell traction forces are generated by actomyosin units within the cell, each of which acts as a force dipole. Therefore, the disks in our model simulations could represent parts of a cell, and their motion represents the dynamics of cell protrusions. The resulting deformation induced by a force dipole in the elastic substrate is given by the strain, $u_{\mathit{ij}}$, which is determined by a force balance in linear elastic theory (SI Appendix, section A), and depends on the material properties of the elastic medium, specifically, the stiffness or Young’s modulus $E$, and the compressibility, given by the Poisson’s ratio $\nu$ (48). The substrate deformation ($u_{\mathit{xx}}$ component of strain) generated by a dipole (oriented along the laboratory x-axis) embedded on the surface of a linear elastic medium is shown in Fig. 1D for two representative values of $\nu$. Here, the blue (red) coloring represents expanded (compressed) regions of the substrate. We note that the extracellular matrix in biological tissue is typically viscoelastic, and over long times, the cell-generated strains may relax. However, our model still applies at short time scales and for linearly elastic synthetic substrates such as polyacrylamide that are routinely used in cell culture experiments (19).

A second contractile force dipole will tend to position itself in and align its axis along the local principal stretch in the medium to reduce the substrate deformation. The resulting interaction potential arises from the minimal coupling of one dipole (denoted by $\beta$) with the medium strain induced by the other (denoted by $\alpha$) and is given by $W^{\alpha \beta }=P_{\mathit{ij}}^{\beta }{u}_{\mathit{ij}}^{\alpha }$ (17). The interaction energy between two dipoles then decays with their separation distance as $W^{\alpha \beta }\sim (P^2/E)·r_{\alpha \beta }^{-3}$. We denote the characteristic elastic interaction energy when the dipoles are separated by only one cell length as, ${\mathcal{E}}_c=P^2/(16E{d}^3)$, where the detailed expression is derived in SI Appendix, section A. This coarse-grained description abstracts out the biophysical details of mechanotransduction but provides a simple physical model for the cell response to deformations in their elastic medium (18).

Representative simulation snapshots (Fig. 1E) of final configurations show that elastic dipolar interactions induce network formation in a stiffness-dependent manner. The Central snapshot corresponds to an optimal substrate stiffness $E^{\ast }$ at which elastic interactions are maximal, while those to the Left (Right) correspond to substrates that are too soft (stiff) for connected network formation. The origin of this optimal stiffness lies in the adaptation of cell contractile forces to their substrate stiffness, as we discuss later.

We expect the multicellular structures resulting from the dipolar cell–cell interactions to depend on three crucial nondimensional combinations of model parameters: the ratio of a characteristic elastic interaction energy ${\mathcal{E}}_c$, to noise – denoted by $A={\mathcal{E}}_c/k_B{T}_{\text{eff}}$ – the effective elastic interaction parameter; the number of cells $N$, equivalently expressed as a cell density or packing fraction, $\varphi =\frac{\pi N{d}^2}{4L^2}$; and Poisson’s ratio, $\nu$, which determines the favorable configurations (both position and orientation) of a pair of dipoles. To show the types of multicellular structures that result from our model elastic interactions, we perform Brownian dynamics simulations (detailed in Materials and Methods) to generate representative snapshots at slices of this $A-\varphi$ parameter space for two values of $\nu$: $0.5$ and $0.1$, shown in Fig. 2 and SI Appendix, Fig. S4, respectively. As packing fraction is increased, networks form more readily. As the effective elastic interaction is increased, cells form into networks characterized by chains, junctions, and rings. This can be thought of naturally as a competition between entropy and energy. At low packing fractions or effective elastic interaction, cells are either in a gas-like state or form local chain segments with many open ends which have high entropy. As packing fraction or effective elastic interaction increases, cells relinquish translational and rotational freedom for more energetically favorable states such as longer chains, junctions, or rings. This is consistent with the cell density–dependent morphologies seen in images from in vitro hydrogel experiments (reproduced from ref. 43) and shown in Fig. 2.

We choose two representative values of $\nu$ in our model simulations because their corresponding strain plots are qualitatively different (37) as seen in Fig. 1D. Briefly, since contractile dipoles prefer to be on stretched regions of the substrate, the low (high) $\nu$ deformation patterns are expected to favor two (four) nearest neighbors. The different values of the Poisson ratio could correspond to synthetic hydrogel substrates and the fibrous extracellular matrix, respectively. While hydrogel substrates are nearly incompressible ($\nu =0.5$), the ECM comprises of networks of fibers which permit remodeling and poroelastic flows leading to reduced material compressibility (e.g., $\nu =0.1$) at long time scales (49).

To characterize the extent of multicellular network formation, we consider the percolation order parameter which quantifies the ability of a connected network to span the available space. Percolation is defined as the probability that, for a steady state realization of the network, there exists a continuous path through it that spans the length of the simulation box. To compute percolation probability, we first identify connected clusters of cells, a process detailed in SI Appendix, section E. A specific network configuration is considered to be percolating if any two cells within the same cluster are separated by a Euclidean distance greater than or equal to the simulation box size. The average values and corresponding errors are then plotted against varying packing fraction $\varphi$ in Fig. 3A and varying effective elastic interaction parameter $A$ in Fig. 3B. Multiple such simulations are then combined into a phase diagram in $A-\varphi$ parameter space in Fig. 3C. The results show that percolation requires both density and interaction strength to be above corresponding threshold values.

To contrast with the dipoles that mutually align through long-range and anisotropic interactions, we consider a control system of “diffusing sticky disks.” These agents just diffuse without any long-range interactions and cease movement upon contact with another agent. We find percolating networks for both interacting elastic dipoles and diffusing sticky disks. However, Fig. 3A shows that model cells which interact as dipoles at long-range require far fewer cells to percolate than their sticky disk counterparts given that the elastic interaction strength is sufficiently greater than noise as shown in Fig. 3B ($A⪆1$ in the case shown where $N=300$). This is because the anisotropic nature of the dipolar interactions promotes end-to-end alignment of cells, creating elongated structures like chains, which can self-assemble into space-spanning networks. We therefore show that network formation requires fewer cells when cells can sense, move, and align in response to the substrate deformations created by other cells. Thus, networks guided by mechanical interactions are more cost efficient than when cells move or spread randomly, forming adhesive contacts upon finding their neighbors.

Much work has been done on characterizing the connectivity percolation transition on various lattice configurations (50). The critical packing fraction can be widely different depending on the lattice geometry and whether the space-spanning clusters comprise sites or bonds (51, 52). The critical packing fraction for site percolation is known to be ${\varphi }_C=0.5$ for an infinitely large triangular lattice (53). In approximate agreement with this, we find that the critical packing fraction for diffusive sticky disks for the current finite system size $L$ is ${\varphi }_C\approx 0.6$. For the dipolar particles, anisotropic interactions shift the percolation transition to ${\varphi }_C\approx 0.2$, similar to those seen in dipolar colloidal assemblies at low reduced temperature (54).

Our observed packing fractions for transition to percolation are specific to the simulation system size, $L$, and differ from the actual critical packing fraction due to finite size effects. How prominent these effects will be depends on the fractal dimension, which provides a measure of how these structures scale with size. Since area scales like $L^2$, but number of particles scales like $L^{d_f}$, where $d_f$ is the fractal dimension, ${\varphi }_C\propto L^{d_f-2}$. Therefore, there exists a regime in which ${\varphi }_C$ will decrease with increasing $L$, as shown by simulations with bigger box sizes (SI Appendix, section F). We present an analysis of the fractal dimension of these networks and corresponding experiments in the next section.

We showed in the previous section that the cells’ ability to form networks is expected to depend on the strength of elastic interactions arising from their mutual deformations of the substrate. To compare with experiments, we now consider how this elastic interaction strength $A$ depends on the substrate stiffness. Experiments show that cells spread and polarize more on substrates of increasing stiffness, such that their traction force saturates to a maximal value $P_0$ at a characteristic substrate stiffness, $E^{\ast }$, that depends on cell type and matrix mechanochemistry. The effective elastic interaction parameter, $A$, can be mapped to substrate stiffness, $E$, by using a model relation predicting the dependence of cell traction force on substrate stiffness (55): $P=P_{0}E/(E+E^{\ast })$. The resulting elastic interaction parameter, $A$, is weak on softer substrates where cell forces are low and also on stiffer substrates, where the deformations are low. It reaches a maximum at the characteristic stiffness $E^{\ast }$ as detailed in SI Appendix, section G. This mapping from effective elastic interaction to substrate stiffness (SI Appendix, Fig. S6) results in a peak in the percolation curves (SI Appendix, Figs. S8 A and C) over an interval of substrate stiffness centered around the optimal stiffness $E^{\ast }$. This interval depends on both cell density and effective temperature representing noisy cell movements. Higher effective temperature and lower cell density reduce both peak height and width. This result is consistent with experiments on EC cultures (Fig. 1A) which show that percolating networks form only in a certain range of substrate stiffness, but these previous works do not demonstrate that network formation is optimal at intermediate substrate stiffness (43, 44).

To test this prediction of our model, we performed 2D cell culture experiments on elastic substrates over a wide range of stiffness values. Human umbilical vascular endothelial cells (HUVECs) were cultured at three different seeding densities ($8\times {10}^3$/${\text{cm}}^2$, $14\times {10}^3$/${\text{cm}}^2$, and $20\times {10}^3$/${\text{cm}}^2$) on fibronectin-coated polyacrylamide substrates of varying stiffness: (200 Pa, 480 Pa, 1 kPa, 2 kPa, 4.5 kPa, and 10 kPa). The substrate preparation protocol, described in Materials and Methods, and stiffness characterization of these substrates follow standard precedents (56). Cells were fluorescently labeled and imaged at regular intervals over the course of 19 h post-seeding. While for the lower seeding density, network formation could be observed at these longer time scales (Fig. 4A, Middle), the higher seeding density led to denser, isotropic clusters, and a resulting loss of network morphology (Fig. 4A, Right). We then considered the images of these denser cultures at 9 h instead of 19 h, where network morphology was still apparent. We also observed that the dense isotropic clusters were more prevalent at higher substrate stiffness due to enhanced cell spreading and possibly proliferation at later times.

To quantitatively obtain the percolation probability for the observed cell clusters, we process the experimental images (Fig. 4A) by emphasizing intercellular connections as described in Materials and Methods under Image Analysis. We then parse the resulting binary images (Fig. 4B) into $N_B=$ 312 subboxes each so as to obtain sufficient statistics from a single experimental image. We next computed the mean percolation probability over all subboxes, $p=\frac{1}{N_B}\sum _{i=1}^{N_B}{p}_i$ and the corresponding standard error of mean. Here, we set $p_i=1$ if the $i$th subbox is “percolating,” i.e., it contains a cluster that spans the subbox, and set $p_i=0$ otherwise. To compare the sparser, heterogeneous experimental configurations with our simulated networks, we normalized these values by the maximum mean percolation probability across all experimental seeding densities and stiffnesses. For practical convenience, we henceforth denote the normalized percolation probability value as $p$.

We find that for both the lowest seeding density sampled at long times (Fig. 4C, Left) and for the highest seeding density sampled at short times (Fig. 4C, Right), the normalized experimental percolation probability exhibits a peak at a stiffness of about 4.5 kPa. Unlike simulations where packing fraction and elastic interaction are independent parameters, the area covered by cells in experiments depends on stiffness because cells spread more on stiffer substrates. This is why we need a range of packing fraction values from simulation to compare with experiment. Like the experimental images, the simulation images were skeletonized to emphasize inter-particle connections (Materials and Methods and Image Analysis). We denote the corresponding packing fraction of skeletonized images by $\tilde{\varphi }$ to distinguish from the packing fraction of simulated disks, $\varphi$. We then plot a family of interpolated simulation curves as a function of substrate stiffness over a range of packing fraction values, $\tilde{\varphi }$, chosen to fit the experimental data in Fig. 4C. These values are close to the range of packing fraction values in experimental images (0.05 to 0.15 for $\frac{8k}{{\text{cm}}^2}$ and 0.1 to 0.2 for $\frac{20k}{{\text{cm}}^2}$). The quantitative agreement of the experimental data with simulation values lends credence to our model that network formation is driven by substrate-mediated elastic interactions and that these are stronger within a range of substrate stiffness values centered around an optimal value, $E^{\ast }$.

We note an important distinction between the predictions of the cell dipole model and the observed cell clusters in experiments. These latter tend to exhibit isotropic dense clusters on stiffer substrates at higher seeding density. We expect this is because cells spread more on stiffer substrates and form direct adhesive contacts with neighbors. Cell spreading and direct cell–cell contact-based interactions are not implemented in our minimal model since we focus on the long-range substrate-mediated dipolar interactions expected to dominate in dilute cultures. At higher densities, endothelial cells are known to form confluent monolayers (57). At intermediate densities, some of these dense isotropic clusters occur alongside networks and elongated structures. Modeling these would require a combination of cell–cell and cell–substrate forces. Although dense isotropic clusters are not seen in our dipole simulations (Fig. 1E), their occurrence in experiment supports our model expectation that the dipolar elastic interaction strength ($A$) becomes smaller on stiffer substrates in relation to the isotropic, cell–cell contact interactions.

While the percolation analysis shown in Fig. 4 validated our model predictions for the substrate stiffness–dependence of network formation, we now seek to predict characteristic morphological traits of the cell clusters. A careful examination of experimental images in Figs. 4 and 5 reveals distinct morphologies of cell clusters, ranging from isolated cells and isotropic clusters to networks and elongated clusters. To obtain a measure of how elongated each cell cluster is, we calculate a “shape parameter,” defined as $s\equiv \frac{R_g^2}{\mathit{Area}}=\frac{1}{N^2}\sum _{k=1}^N{({\mathbf{r}}_k-{\mathbf{r}}_{\mathit{CM}})}^2$, for each unique cluster as described in Materials and Methods under Image Analysis. Here, $R_g$ represents the radius of gyration of the cluster, which is defined about its center-of-mass ${\mathbf{r}}_{\mathit{CM}}$, and $N$ is the number of occupied pixels in each cluster. The normalization by cluster area ensures that we control for cluster size variations between different experiments. Lower values of this shape parameter correspond to isotropic shapes, the theoretical lower bound being $\frac{1}{2\pi }$ for a solid circular disk. Conversely, a higher shape parameter corresponds to more elongated clusters. To compare with simulation, we scale the shape parameter of each cluster by their global mean across all identified clusters at the different seeding densities and stiffnesses. Henceforth, the scaled values of the shape parameter are denoted by $s$.

To classify the dominant morphological feature in each image, we constructed a composite order parameter combining the global information of connectivity percolation with the local cluster-scale morphological characteristics captured by the shape factor. The order parameter is defined such that clusters with normalized percolation probability above a threshold value ($p_T=0.7$) are considered “networks.” We choose this value of $p_T$ to pick out experimental images where a few of the largest clusters contain more than 20% of the total filled area (SI Appendix, Fig. S20). If $p<p_T$, implying that there are no dominant space-spanning clusters, we classify clusters into “isolated” or “chains,” depending on whether $s$ is less or greater than a threshold value $s_T=0.95$. This value of $s_T$ is chosen to correspond to simulations with two aligned dipoles, giving an elongated morphology that this parameter is designed to capture. The order parameter which accomplishes the above classification is given by $\text{OP}\equiv \mathrm{\Theta }(p-p_T)p+\mathrm{\Theta }(p_T-p)(0.25+0.5(s-s_T))$, where $\mathrm{\Theta }(x)$ is the Heaviside step function and the numerical factors give an $0<OP<1$ for the specific values of $s_T$ and $p_T$, justified above. The differences are captured by ranges of values of the order parameter: $OP>0.7$, 0.25 to 0.7, and $<0.25$ correspond to percolating “networks,” elongated but disconnected “chains,” and isotropic “isolated” clusters, respectively.

We compute this order parameter for interpolated simulation data and for experimental data 9 h post seeding, so that cell proliferation effects are minimal. Experimental data once again reveal a nonmonotonicity in network formation at the high densities in Fig. 5, Bottom Left. We construct a phase diagram of this order parameter in $\tilde{\varphi }-E$ space (Fig. 5, Center), where the color map corresponds to simulation data. Overlaid on this phase diagram are discrete markers representing experimental data, which have been classified into the three distinct regimes according to their measured order parameter values. The data comprise three initial seeding densities and six stiffness values giving a total of eighteen data points. Their distribution clearly shows a correlation between cell area coverage and substrate stiffness. This is due to cells spreading more on stiffer substrates, readily seen through the lack of low packing fraction data at higher stiffness.

Phase boundaries are drawn as dashed lines that delineate the distinct regions of the simulation order parameter values. Of the eighteen experimental data points, only two lie outside of the corresponding predicted regions. Both of these are at low substrate stiffness and intermediate packing fraction. These are classified as “chains,” but lie in the “isolated” part of the predicted phase diagram. We expect that at these intermediate densities, cells can spread and touch each other to form elongated structures even if the elastic dipolar interactions are small. Since our model does not include such spreading effects, this discrepancy is not surprising. Overall, the model phase diagram closely predicts the experimentally observed multicellular structures.

We now further evaluate the morphological similarity of the networks from our simulated dipoles and our experimental cell culture by calculating the fractal dimension. For the “sticky disks,” we find a fractal dimension of $d_f=1.81$, whereas for the dipoles, we find fractal dimensions of $d_f=1.698\pm 0.004$ and $d_f=1.711\pm 0.003$ for $\nu =0.1$ and $\nu =0.5$, respectively. We find a similar fractal dimension for our experimental HUVEC culture in the network regime on a substrate of stiffness $E=4.5\text{kPa}$, $d_f=1.722$. Interestingly, simulated networks on substrates of $\nu =0.1$ and $\nu =0.5$ are statistically distinguishable, with the experimental fractal dimension showing reasonable agreement with the $\nu =0.5$ simulated dipole case. This is in accordance with the approximately incompressible nature of hydrogel substrates. The proximity of the fractal dimensions of the simulated dipoles to that of experimental cell networks, in relation to the sticky disks, indicates that cells utilize a more complex strategy to self-assemble than simply random movement followed by cell–cell adhesion. The elastic dipolar interactions are thus a plausible strategy allowing the self-assembly of biologically desirable, space-spanning, and cost-effective networks.

We now focus on simulated networks (such as in Fig. 3) to thoroughly characterize their two predominant network morphological constituents—branches and rings. We relate the resulting structural metrics to the transport function of biological networks. We highlight qualitative differences in morphology of simulated networks between elastic substrates with high and low Poisson’s ratio values, which may motivate future experimental investigation. Fig. 6A shows average branch length for $N=300$ ($\varphi =0.33$) cells as a function of effective elastic interaction ($A$). The average branch length for the higher $\nu$ case remains roughly constant and low at about two cell lengths. The lower $\nu$ case exhibits a peak in average branch length at the percolation threshold ($A=1$) before decreasing and saturating at high $A$ values. The distribution of branch lengths (Fig. 6B) shows that while $\nu =0.5$ is sharply peaked at $d$, $\nu =0.1$ exhibits branches greater than 18 $d$ and shows a greater relative count in the 3 to 10 $d$ range.

These results suggest that at higher values of $\nu$, network morphology is more resilient to noise, and the branch lengths are not as easily tunable, The greater variability in branch lengths leads to longer branches in the lower $\nu =0.1$ case, which then requires (for $A\ge 5$) fewer cells to percolate than at $\nu =0.5$. This is seen by the difference of the curves at the shoulder of the percolation transition in Fig. 3A and SI Appendix, Fig. S10. The greater resilience of the network at higher substrate $\nu$ leads to percolation at smaller $A$ than its low $\nu$ counterpart (Fig. 3B and SI Appendix, Fig. S10). In SI Appendix, section I, we construct a detailed map of the percolation transition in the $A-\varphi$ parameter space, to show how $\nu =0.1$ requires fewer cells to percolate for a range of $A$ values, while $\nu =0.5$ can percolate at lower values of $A$. This suggests that the two regimes of substrate compressibility optimize two different measures of cost of network building: one, the number of cells, and the other, the strength of cell contractility.

Fig. 6C shows a cumulative distribution of ring area for our networks at two crucial regions in our parameter space—those at which the networks are well above the percolation transition (solid lines), or just above it (dashed lines). Similar to the branch length distribution, the networks at higher $\nu$ form many small rings and few large rings, while the lower $\nu$ case shows a broader distribution of ring sizes. The tendency of the $\nu =0.5$ configurations to form numerous smaller rings leads to a marginally less efficient area coverage than the low $\nu$ case, which forms longer branches and fewer small rings (SI Appendix, Fig. S11). These results are also consistent with the fractal dimensions we obtained earlier, with $d_f$ being slightly higher for the $\nu =0.5$ than the $0.1$ cases, indicating more compact structures for the former. These topological features also give rise to distinct coordination numbers for the two compressibility regimes (SI Appendix, Fig. S18). Interestingly, the coordination number on the lower Poisson ratio substrate resembles those near the rigidity percolation of elastic fiber networks (58, 59)—indicating connectivity percolation as a precursor to mechanical rigidity which is relevant to both tissue development and disease (60).

To examine the robustness of our model networks to damage, a biologically significant property, we measure the largest remaining cluster size as a function of the fraction of network bonds removed (61) (Fig. 6D). We find that whether well above or just at the percolation threshold, the networks at higher $\nu$ retain cluster size well as bonds are removed. Networks at lower $\nu$ well above the percolation transition lose largest cluster size at the same rate as their higher $\nu$ counterpart. At the shoulder of percolation, however, networks at low $\nu$ lose largest cluster size and fall apart much more rapidly than any of the other networks. This is the same parameter regime at which networks at low $\nu$ exhibit a peak in branch length. By forming long branches, ring structure formation is sacrificed. Thus, we find that the prime factor for robust networks is the tendency to form rings which provide degeneracy to paths between any two nodes in the network—a result consistent with network structure optimization models (62). In summary, at lower $\nu$, networks tend to form longer and more broadly distributed branches which promote efficiency with respect to the filling and spanning of space at the cost of being susceptible to damage, while at higher $\nu$, networks are predominantly composed of small rings, which provide robustness to the networks at the cost of transport efficiency.

We model the ubiquitous traction force pattern of a polarized cell as a single, anisotropic force dipole. The dipole magnitude is the cell force times the distance along the long axis of the cell, $P=Fa$. Since the contractile cytoskeletal machinery (e.g., actomyosin stress fibers) of the cell is typically aligned along this axis, this is also usually the principal direction of stress exerted by the cell and is henceforth called the “dipole axis.” Such a force dipole induces a strain in the substrate, which is modeled as an infinitely thick, linear, isotropic elastic medium.

By considering two dipoles ${\mathbf{P}}^{\alpha }$ and ${\mathbf{P}}^{\beta }$, we show in SI Appendix, section A that the work done by a dipole $\beta$ in deforming the elastic medium in the presence of the strain created by the other dipole $\alpha$ is given by (37): $W^{\alpha \beta }=P_{\mathit{il}}^{\beta }{u}_{\mathit{il}}^{\alpha }({\mathbf{r}}^{\beta })$, where the strain can be written in terms of ${\mathbf{P}}^{\alpha }$ and second derivatives of an elastic Green’s function as $u_{\mathit{il}}^{\alpha }({\mathbf{r}}^{\beta })={\partial }_l{\partial }_k{G}_{\mathit{ij}}({\mathbf{r}}^{\beta }-{\mathbf{r}}^{\alpha })P_{\mathit{jk}}^{\alpha }$. This minimal coupling between dipolar stress and medium strain represents the mechanical interaction energy between dipoles. Typical substrate strains are shown in Fig. 1D, where the blue (red) coloring represents expanded (compressed) regions. A second or test dipole present in these regions would tend to align its contractile axis along the principal stretch direction of the substrate. In the expanded (blue) regions, the test dipole is aligned with and attracted toward the central dipole, whereas in the compressed (red) regions, a test dipole is aligned orthogonal to and repelled away from the central dipole. The orientational dependence of the strain field is changed by the Poisson’s ratio or compressibility of the substrate (18).

Our computational “many-cell” model considers cells as discrete agents ($N$ agents in a $L\times L$ box with periodic boundary conditions) which move and orient randomly, but that also interact with one another through long-range elastic interactions via a force dipole strain field coupling and a short-range repulsive spring. Fig. 1E shows our simulation setup and the main ingredients of the model. We ignore details of the cell shape and subcellular structures in this minimal model and instead consider the cells as disk-shaped agents endowed with contractile, elastic dipoles. This simplifying assumption implies that we do not consider changes in the shape and size of individual cells that occur as a result of cell–substrate feedback when substrate stiffness is varied but instead focus on the multicellular structures at longer length scales.

We now consider the translational and orientational dynamics of a collection of model cells. These interact with each other through short-range, steric and long-range, substrate-mediated, elastic interactions, and undergo diffusive motion. The overdamped Langevin dynamics governing the position of a cell labeled $\alpha$ is

where $D_{\text{T}}$ is the effective translational diffusivity quantifying the random motion of an isolated moving cell, with ${\mathit{\eta }}_{\text{T}}$ as a random white noise term whose components satisfy $⟨{\eta }_{i,\text{T}}(t){\eta }_{j,\text{T}}(t^{\prime })⟩=\delta (t-t^{\prime }){\delta }_{\mathit{ij}}$. Typical adherent cells do not move very persistently, and at time scales much longer than their persistence time, their motion is random and has been shown to be well characterized by a diffusion constant (69). We thus neglect the directed self-propulsion term typically included for active particles from the dynamics. The mobility ${\mu }_T$ in Eq. 1 is inversely related to the effective friction from the medium that the moving cell experiences at its adhesive contacts with the substrate. Similarly, the orientational dynamics of the cell denoted by $\alpha$ is given by

where ${\widehat{\mathbf{n}}}_{\alpha }$ is the unit vector along the dipole axis of the cell $\alpha$ and $D_{\text{R}}$ is the effective rotational diffusivity quantifying the random reorientations of an isolated moving cell. Cells encounter various forms of internal stochastic effects including internal cytoskeletal rearrangements producing membrane morphological fluctuations, substrate surface binding fluctuations, and fluctuations in myosin motor forces, which are all absorbed into a coarse-grained effective temperature, $T_{\text{eff}}$, in our model. Single cell and cell cluster experiments have shown this effective temperature to be on the order of ${10}^{-15}$ to ${10}^{-14}$ J (70). Though the rotational and translational diffusion are in principle independent, we will here assume them to correspond to the same underlying processes and therefore the same effective temperature, $k_B{T}_{\text{eff}}=D_T/{\mu }_T=D_R/{\mu }_R$. We also show some exceptions to this assumption in SI Appendix, section L, which all robustly form networks.

The pairwise cell–cell interaction potential $W_{\alpha \beta }$ between cells labeled $\alpha$ and $\beta$ consists of the long-range elastic interaction arising through their mutual deformation of the substrate (SI Appendix, section A), and a short-range steric interaction between two cells in contact, and is given by,

where $f$ is a function of Poisson’s ratio—shown in SI Appendix, section A, ${\theta }_{\alpha }$, and ${\theta }_{\beta }$ where $cos{\theta }_{\alpha }={\widehat{\mathit{n}}}_{\alpha }·{\widehat{\mathit{r}}}_{\alpha \beta }$ and $cos{\theta }_{\beta }={\widehat{\mathit{n}}}_{\beta }·{\widehat{\mathit{r}}}_{\alpha \beta }$ are the orientations of cell $\alpha$ and cell $\beta$ with respect to their separation vector, ${\mathbf{r}}_{\alpha \beta }={\mathbf{r}}_{\beta }-{\mathbf{r}}_{\alpha }$ connecting the centers of the two model cell dipoles, respectively. Note that while the elastic potential is in principle long range, it decays strongly as a $1/r^3$ power law; we cut this pairwise interaction off at $r_{\alpha \beta }>3d$ in our simulations, since the substrate strain induced by one cell is unlikely to be detected by a cell few cell lengths away (24).

The above equations are nondimensionalized by a suitable choice of length, time, and energy scales. By choosing the length scale to be the cell diameter $d$, the time scale to be an elastic time, $t_c=\frac{16E{d}^5}{P^2{\mu }_{\text{T}}}$, and a characteristic elastic interaction as the energy scale, ${\mathcal{E}}_c=\frac{P^2}{16E{d}^3}$, the dynamical equations reduce to (Appendix B),

for the translational motion, while the rotational equation of motion can be written as

where the starred variables indicate nondimensionalized quantities and we have assumed ${\mu }_{\text{R}}{d}^2={\mu }_{\text{T}}$ and $D_{\text{R}}{d}^2=D_{\text{T}}$, although the latter is not required for a system that is out of equilibrium. The nondimensionalized pairwise interaction potential in Eq. 3 is here given by $W_{\alpha \beta }^{\ast }=\frac{1}{2}{k}^{\ast }{(1-r^{\ast })}^2\mathrm{\Theta }(1-r^{\ast })-\frac{16f}{r^{\ast 3}}\mathrm{\Theta }(r^{\ast }-1)$, where $k^{\ast }=k{d}^2/{\mathcal{E}}_c$. We introduce an effective elastic interaction parameter quantifying the elastic interaction strength relative to that of intrinsic noise in the cell motion,

where the noisy cell movements correspond to an effective temperature, $k_B{T}_{\text{eff}}\equiv D_T/{\mu }_T$. This explicitly shows that $A$ is a measure of the characteristic elastic interaction energy scale relative to the magnitude of cell stochasticity described by an effective temperature.

In experiments, the value of the effective interaction parameter $A$ will depend on cell contractility, the stiffness of the elastic substrate, and the diffusivity that originates from the motility of single cells. Importantly, cells adapt their contractile forces to the stiffness of the underlying substrate. Measurements (71) and models (55) of the dependence of cell force on substrate stiffness suggest that the magnitude of the force dipole can be written as $P(E)=P_0/(1+E/E^{\ast })$, where the characteristic substrate stiffness for a given cell at which the cell traction forces saturate to their maximal value is denoted by $E^{\ast }$. This dependence when inserted into the definition of the effective elastic interaction parameter, $A$, in Eq. 6, leads to $A$ being a peaked function of $E$. Since stiffer substrates are harder to deform and cells on softer substrates do not generate enough traction, substrate deformations and therefore elastic interactions are maximal at an intermediate optimal stiffness value ($E=E^{\ast }$) (Table 1).

To identify a plausible range for the values of $A$ consistent with cell culture experiments, we note that the typical values for the force dipole for contractile cells adhered to elastic substrates is $P_0=Fd\sim {10}^{-12}$ to ${10}^{-11}$ J (31, 38). This corresponds to measured traction forces of $F\sim$ 10 to 100 nN with a distance of $\sim$50 $\mathrm{\mu }$m separating the adhesion sites at which the forces act on the substrate (29, 69), which is also the typical size of the cell along its long axis. For a typical substrate stiffness of $E\sim 1$ kPa characteristic of EC network formation (43, 44), we therefore estimate an elastic dipole energy of ${\mathcal{E}}_c=\frac{P^2}{16E{d}^3}=\frac{F^2}{16Ed}\sim {10}^{-15}$ J, similar to measured values for cell contractile energy stored in the elastic substrate (72). Since adherent cells crawl by exerting forces at the focal adhesions at which forces are transmitted to the substrate, the net mobility that determines cell translation, ${\mu }_T$, can be estimated from the friction force at these adhesion sites. From the observation that the focal adhesions with surface area of $10$ $\mathrm{\mu }$m² reorient with speeds of $\mathrm{\mu }$m/min in the direction of an external, applied stress of kPa (73), we can estimate the mobility coefficient (inverse of friction coefficient) to be ${\mu }_T\sim 0.1\mathrm{\mu }$m/min pN${}^{-1}$. The effective diffusivity characterizing single endothelial cell movements was measured to be $\sim$10 $\mathrm{\mu }{\text{m}}^2/$min (23, 43, 74). Together, these give an estimate for the effective temperature: $k_B{T}_{\text{eff}}=D_T/{\mu }_T\sim {10}^{-16}$ J $\sim {10}^4$$k_{B}T$. For substrate stiffness $E\sim E^{\ast }1$ kPa, we thus estimate the ratio of elastic energy to noise to be $A={\mathcal{E}}_c/k_B{T}_{\text{eff}}\sim 10$.

In experiments, the substrate stiffness can be tuned over a wide range. In particular, Califano et al. tested the formation of EC networks on substrates whose rigidity was varied from 100 Pa to 10 kPa (44). This, in our estimate, corresponds to an interaction parameter $A\sim 1-100$, with $A=0.1$ corresponding to high noise or nonoptimal values of substrate stiffness (too soft or too stiff). Similarly, we can estimate the characteristic timescale as $t_c=\frac{d^2}{{\mathcal{E}}_c{\mu }_{\text{T}}}\sim {10}^2$ min. This timescale of hours is consistent with that required for the formation of cellular structures in experiments (44).

The following processing is done in order to directly compare simulation predictions to experimental results (Figs. 4 and 5) and to obtain network metrics for simulations (Fig. 6 and SI Appendix, Figs. S9, S11, S12, and S17). All image analysis used in this work was carried out using the open-source software ImageJ (75). Raw grayscale experimental images are imported into ImageJ. “Enhance Contrast” command is run with “saturated pixels” widget set to 2. We then “Despeckle” the image and “Enhance Contrast” once more before running “Subtract Background” with a rolling ball radius of 50 pixels. We “Gaussian Blur” with a sigma of 10 pixels. We then threshold, keeping intensities 20 and above. This is then converted into a mask, skeletonized, and dilated four times so as to preserve the raw filling fraction (Fig. 4B).

For simulated networks like those shown in Fig. 2, we replace the isotropic disk markers with “pill-box” shaped markers (as seen in SI Appendix, Fig. S11) which are elongated along the dipole axis of each cell to guide the subsequent skeletonization. Using ImageJ, we first apply “Gaussian Blur” with a sigma of 2 pixels; then, we threshold keeping intensities 150 and above and then convert into a mask and skeletonize. Finally, we dilate the skeleton four times so that small-scale features of assembly like compact rings are preserved, while washing out the shape of the individual disks. At this point, both experimental and simulated images are dilated skeletons. The packing fraction of the dilated skeleton representations of simulations ($\tilde{\varphi }$) are smaller than their respective nonoverlapping disk packing fraction ($\varphi$) by a factor, $\tilde{\varphi }\approx 0.75\varphi$. To compute fractal dimensions, we follow the aforementioned image processing with the additional step of dilating experimental skeletons so as to have roughly the same packing fraction as simulations We then use ImageJ’s “Fractal box count” tool with the default pixel array.

To identify unique clusters in both experimental and simulation images, dilated skeleton images are imported into a custom python program. This program assigns a cluster label to the first nonzero pixel and then does recurrent loops assigning neighboring pixels to the same cluster label until a pixel is identified that does not neighbor any of the pixels with this cluster label. The cluster label is incremented and the process repeats until every nonzero pixel is assigned a cluster label.

P.S.N., J.E.Z.A., K.E.M., A.G., and K.D. designed research; P.S.N., J.E.Z.A., F.G., A.G., and K.D. performed research; J.E.Z.A., F.G., and K.E.M. contributed new reagents/analytic tools; P.S.N. and J.E.Z.A. analyzed data; and P.S.N., K.E.M., A.G., and K.D. wrote the paper.

The authors declare no competing interest.