Microstructural properties of extracellular matrix (ECM) promote cell and tissue homeostasis as well as contribute to the formation and progression of disease. In order to understand how microstructural properties influence the mechanical properties and traction force-induced remodeling of ECM, we developed an agent-based model that incorporates repetitively applied traction force within a discrete fiber network. An important difference between our model and similar finite element models is that by implementing more biologically realistic dynamic traction, we can explore a greater range of matrix remodeling. Here, we validated our model by reproducing qualitative trends observed in three sets of experimental data reported by others: tensile and shear testing of cell-free collagen gels, collagen remodeling around a single isolated cell, and collagen remodeling between pairs of cells. In response to tensile and shear strain, simulated acellular networks with straight fibrils exhibited biphasic stress–strain curves indicative of strain-stiffening. When fibril curvature was introduced, stress–strain curves shifted to the right, delaying the onset of strain-stiffening. Our data support the notion that strain-stiffening might occur as individual fibrils successively align along the axis of strain and become engaged in tension. In simulations with a single, contractile cell, peak collagen displacement occurred closest to the cell and decreased with increasing distance. In simulations with two cells, compaction of collagen between cells appeared inversely related to the initial distance between cells. These results for cell-populated collagen networks match in vitro findings. A demonstrable benefit of modeling is that it allows for further analysis not feasible with experimentation. Within two-cell simulations, strain energy within the collagen network measured from the final state was relatively uniform around the outer surface of cells separated by 250 *μ*m, but became increasingly nonuniform as the distance between cells decreased. For cells separated by 75 and 100 *μ*m, strain energy peaked in the direction toward the other cell in the region in which fibrils become highly aligned and reached a minimum adjacent to this region, not on the opposite side of the cell as might be expected. This pattern of strain energy was partly attributable to the pattern of collagen compaction, but was still present when mapping strain energy divided by collagen density. Findings like these are of interest because fibril alignment, density, and strain energy may each contribute to contact guidance during tissue morphogenesis.

## Introduction

The extracellular matrix (ECM) provides rich biochemical and biomechanical cues that regulate many cell processes such as adhesion, migration [1–3], proliferation [4], differentiation [5,6], apoptosis [7], and intercellular signaling [8]. In vivo, cells continuously maintain their local ECM through a balance of matrix synthesis, chemical modification, degradation, and the application of traction forces [9,10]. It is believed that a breakdown of this balance and subsequent changes to the mechanical properties of the ECM could disrupt homeostasis and promote pathogenesis [11] such as cancer metastasis [4,12] and fibrosis [13]. Given the importance of the mechanical properties of the ECM for promoting and maintaining a desirable cell phenotype, a greater understanding of how cells locally remodel the ECM and how the microstructural properties of the ECM influence this remodeling has implications for the treatment and prevention of many different diseases.

Reconstituted type I collagen hydrogels are commonly used scaffolds for studying ECM remodeling. Collagen is a natural biomaterial possessing a fibrous microstructure that is biophysically and biochemically familiar to cultured cells. Additionally, the fibrous nature of ECM proteins like collagen is critical both for remodeling and stress transmission [14–19]. In vitro models using collagen gels have greatly contributed to our understanding of ECM remodeling. Unfortunately, reconstituted type I collagen and other fibrous matrices have inherent limitations. This is because it is very challenging to alter specific properties independent of all others. For example, collagen stiffness (i.e., Young's modulus or storage modulus) is known to influence cell behavior and is often changed by varying collagen concentration [20]. Collagen concentration, however, can also influence fibril number, fibril diameter, porosity, diffusion rates, and the elastic modulus of individual fibrils [21–23]. Thus, in this example, it is not clear whether variations in cell response are due to stiffness or the influence of other variables that also change with concentration. Of further consideration, bulk mechanical properties (e.g., storage or loss moduli) do not accurately represent the complex biophysical and biochemical structure of the microenvironment that a cell experiences or allow for the accurate prediction of remodeling that will occur in response to cell traction force. These reasons motivate the study of how individual microstructural parameters (e.g., fibril modulus, diameter, interconnectivity, and availability of integrin binding sites) influence cell–matrix interaction, remodeling, and cell migration.

Computational modeling is complementary to the reconstituted type I collagen in vitro model and may provide insights otherwise unattainable with experimental systems. One computational approach that possesses fibril-level detail is the structural model in which the overall behavior of the network results directly from the summed responses of individually modeled collagen fibrils and cross-links [15,24–29]. This approach provides the ability to probe the collagen microstructure at a fine level of spatial and temporal detail and it allows for a rigorous analysis of the sensitivity to individual system properties. Using finite element-based structural models, researchers have demonstrated how collagen alignment contributes to the macroscopic behavior of collagen gels [30] and that strain-stiffening may result from geometric realignment of the network as opposed to strain-stiffening of individual fibrils [24,25]. Finite element models have also been used to predict collagen fibril orientation in collagen constructs in the presence/absence of cell traction force [24,26]. In addition, this modeling technique has been used to estimate the radius within which fibrils align around a contracting cell [27]. Working with a two-dimensional (2D) random network, Abhilash et al. expanded upon this idea and explored how cell contractility and cell aspect ratio contribute to a “zone of mechanical influence.” They also explored how separation distance between a pair of cells influences long distance, mechanically mediated, intercellular communication along bundles of aligned fibers [28]. Humphries et al. demonstrated that the initial fiber network geometry can influence the extent of fiber alignment, substrate displacement, and stress propagation [31]. Despite these insights, finite element analysis may have limitations. In the examples that included cell traction, traction was simulated by shortening fibrils within a given region [26], shortening “pseudopodia” attached to the collagen network [27], or contracting an elliptical void representing a cell by a percentage of its area [28,31]. In each case, traction was implemented as a single pull of the cell on the network. This is in stark contrast to the repetitive “hand-over-hand” interactions seen experimentally for cells within collagen gels; cells extend pseudopodia that each may engage a collagen fiber, pull the fiber, release, and repeat [32]. One important impact of modeling cell traction as a “single pull” opposed to more realistic “many pulls” is the magnitude of cell-induced network deformation the model can produce. In each of the finite element models mentioned previously, the change in area or volume of the fibrous network due to cell traction was modest, and much less than the more than 96% decrease in area observed experimentally for cell populated collagen gels [33].

In order to include dynamic cell traction and fibril-level detail, we developed an agent-based model of cell–matrix interactions. This approach has been used to model other dynamic biological processes including angiogenesis in a fibrous ECM as well as cell migration on a surface possessing a stiffness gradient and in a collagen network [34–36]. In agent-based models, agents are represented by independent nodes within a simulated environment. The behavior of these agents and their interactions with other agents are governed by user-defined rules. From a global perspective, the interactions among agents following these simple rules results in complex systemic patterns, often described as “emergent behavior,” which were not intentionally programmed in the model [37]. We have previously shown that our model recapitulates several qualitative aspects of collagen network remodeling such as global compaction and the formation of matrical tracts-relatively dense and aligned fibrils between pairs of cells [15,38]. While the ability of the agent-based model to recapitulate some of the behaviors commonly seen in cell-populated gels is promising, there is a need to better validate this model. In contrast to finite element models, in which the underlying computational approaches as well as their potential applications and limitations have been extensively studied, agent-based mechanical models are relatively poorly characterized and few models have been validated against experimental data.

As a first step toward this validation, we previously reported that our agent-based model accurately predicts single fibril mechanics [15]. Here, we (1) validate our model by comparing simulations results to experimental data reported by others and (2) use the model to gain insight into these experimental systems. To do so, we have selected three sets of experimental data: tensile and shear testing of cell-free collagen gels [24], remodeling around a single isolated cell [39], and remodeling between pairs of cells [40]. It is important to point out that there are differences between the experimental systems used to collect these three data sets including differences in the conditions during collagen polymerization, which are known to profoundly affect collagen microarchitecture and mechanics [23,41], as well as the use of different cell types, which generate different amounts of traction force [42]. In addition, matrix stiffness and collagen concentration have been shown to influence the magnitude of traction force exerted by a single cell type [42]. For these reasons, we are not interested in reproducing specific numerical results, but rather the general patterns of remodeling contained within these data. Therefore, we instead judge validation by the qualitative reproduction of trends in the experimental data.

## Materials and Methods

### General Model Components and Notable Changes to the Model.

For a more complete description of the model presented here, we refer the reader to our prior publications [15,38]. Briefly, the simulations performed here were constructed in NetLogo [43].^{2} All simulations presented here are 2D representations of a 1 *μ*m thick slice through a collagen gel or a cell populated collagen gel; thus, 2D representations of three-dimensional (3D) biological systems and phenomena. In this model, both the cells and collagen fibrils were broken into multiple agents that interacted based on physically realistic rules. Cells were each composed of 109 agents (37 nodes and 72 links). Of the 37 nodes, 30 compose a cell membrane 20 *μ*m in diameter and seven were used to create a central nucleus 10 *μ*m in diameter. These nodes were interconnected with links acting as Hookean springs as shown in Fig. 1 to allow forces to be transmitted intracellularly [38]. The mass for each node was estimated to be 8.484 × 10^{−15 }kg determined by dividing the mass of the cell represented (with a diameter of 20 *μ*m, a thickness of 1 *μ*m, and assuming a density of 1 g/cm^{3}) by the number of nodes (i.e., 37).

Individual fibrils were represented as series of nodes initially 5 *μ*m apart. As described previously, rules were written for the interactions between neighboring nodes so that individual fibrils would behave similar to elastic rods [15]. In stretching and compression fibrils behaved according to Hooke's Law (Eq. (S2), which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection). Bending of fibrils was modeled using torsional springs (Eq. (S3), which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection). Using this approach, we created fibrils that closely matched theoretically predicted deformation for Euler–Bernoulli beams across a wide range of fibril lengths and loading conditions [15]. To form a fiber network, fibrils were added until the target collagen density was achieved. When two fibrils crossed, a cross-link was added, which was modeled as a torsional spring as done by Sander and our previous work [27,38]. Cross-linking reduced the length of links involved to a value between 0 and 5 *μ*m. Surrounding nodes were then adjusted so that link lengths ranged between 2.5 and 5 *μ*m.

Nodes composing the cell membrane possessed the ability to extend a pseudopodium-like link and bind randomly to a node on the collagen network between 3.25 and 6.50 *μ*m away. Each pseudopodium was modeled to pull with a constant force of ∼14 nN. Contraction pulled the fibril toward the cell and, to a very small degree, the cell toward the fibril. After a pseudopodium pulled the node on the network within 1.5 *μ*m from the membrane node, the pseudopodium would release and the membrane node could then extend a new pseudopodium. This hand-over-hand process was repeated indefinitely for all 30 nodes on the cell membrane. If all 30 pseudopodia were pulling simultaneously, a simulated cell could generate a peak traction force of 424 nN, within the wide range of values for total cell traction force reported in the literature (0.1–1500 nN) [44].

Major changes to the model relative to previously published version are summarized in the Supplemental Material, which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection, and include changes in model parameters, an improved cross-linking algorithm. Supplemental Material, which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection, also contains a more complete discussion of the mechanical rules.

For the following models, the dimensions of the region simulated were chosen with the experimental system used for validation in mind. In each, two competing constraints were considered. First, the simulated network needed to be large enough so that deformation would not occur at the edge so that simulation results would be independent of whether the edges were free-floating or constrained. Second, simulation areas needed to be as small as possible to relieve the computational burden imposed by an excess number of agents.

### Three Types of Models

#### No Cell Model.

Randomly distributed fibril networks measured 105 × 105 *μ*m and represent 2.0, 3.0, and 4.0 mg/mL collagen (*n* = 10). Fibrils that do not contribute to the mechanical properties of the networks (i.e., fibrils having dead ends that do not terminate at a fixed edge of the network) were removed to increase computational efficiency (Figs. 2(a) and 2(b)). Differences in collagen concentration were modeled as differences in the number of fibrils. 3.0 mg/mL networks were constructed by adding a number of fibrils equivalent to 1.0 mg/mL to each 2.0 mg/mL network (Fig. 2(c)). Similarly, 4.0 mg/mL networks were constructed by adding fibrils to each 3.0 mg/mL network. In most simulations, fibrils were modeled as initially being perfectly straight. In select simulations, fibril curvature was added to these networks by moving nodes between cross-links perpendicular to the orientation of fibrils (i.e., 0.5, 1.0, 1.5, 2.0, or 2.5 *μ*m) while leaving the nodes at cross-links fixed (Fig. 2(d)). This lengthened fibrils and did not introduce internal stress. Since a node could be moved in one of two opposite directions perpendicular to the fibril and moving nodes in random directions would introduce sharp kinks and lengthen fibrils unpredictably, we arbitrarily chose to move nodes in the downward direction (i.e., between 90 deg and 270 deg) as evident in Fig. 2(d). Additional nodes were added where necessary so that links did not exceed 5 *μ*m.

*ε*), the nodes at the left and right edges of the network were displaced outward in the

*X*direction while allowing their

*Y*coordinates to move freely. To simulate shear strain (

*γ*), the nodes on the top and bottom edges were displaced in opposite directions by forcing both their

*X*and

*Y*coordinates. Tensile strain, or shear strain, was varied incrementally between 0 and 1.05. These values of strain, as with all others in this manuscript, are measured as engineering strain, or engineering shear strain, so that these results can be compared to those from Stein et al. [24]. In these simulations, the stress created by each imposed strain was provided 50,000 iterations to propagate so that the network would approach mechanical equilibrium. Afterward, forces were measured at each node along the edges (i.e., left and right for tensile strain, top and bottom for shear strain). A value of 50,000 iterations was chosen to balance the competing concerns of accuracy and computational demand. This rationale is described in detail in the Supplemental Material, which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection. Forces along each edge were summed and divided by the network cross-sectional area to calculate the stress along that edge. Stresses along the left and right edges were averaged to obtain the average stress (

*σ*). Stresses along the top and bottom edges were averaged to obtain the average shear stress (

*τ*). Incremental Young's moduli and incremental shear moduli were calculated using the stresses and strains at two consecutive measurements (

*σ*

_{2}−

*σ*

_{1}/

*ε*

_{2}−

*ε*

_{1}or

*τ*

_{2}−

*τ*

_{1}/

*γ*

_{2}−

*γ*

_{1}, respectively). Strength of fibril alignment was quantified for the entire network by calculating the anisotropy index (AI). First, the orientation (

*θ*) and length (

_{i}*l*) of each fibril segment along with a reference angle (0 deg ≤

_{i}*θ*< 180 deg) were substituted into Eq. (1) (NF, number of fibril segments).

_{r}*θ*was varied to find the value that maximizes the

_{r}*x*-component of the orientation tensor (Ω

*) [26,45]. That particular reference angle represents the principal angle about which fibril segments are most aligned. The anisotropy index was calculated by substituting the maximum value of Ω*

_{xx}*from Eq. (1) into Eq. (2). AI ranges between 0 and 1 with 0 corresponding to an isotropic distribution and 1 corresponding to perfect alignment among all fibril segments*

_{xx}#### One-Cell Model.

Circular, unconstrained collagen networks 1000 *μ*m in diameter were generated at 2.0 mg/mL (*n* = 12). At initial conditions, a single cell was placed at the center, or origin, of this network. The locations of fibril segments were recorded for every 36,000 iterations throughout simulations. Average displacement was calculated for fibril segments that resided in concentric rings around the cell at the beginning of the simulation (i.e., with a radial distance of 100 ± 3 *μ*m, 133 ± 3 *μ*m, 167 ± 3 *μ*m, etc., chosen to match experimental data from Mohammadi et al. [46]).

#### Two-Cell Model.

In one set of experiments, collagen was modeled using a fibril diameter of 60 nm. Randomly distributed fibril networks measuring 805 × 605 *μ*m were generated at three collagen densities (1.5, 2.0, and 3.0 mg/mL). Differences in concentration were modeled as differences in the number of fibrils. Two cells were centered on each network and separated by 75, 100, 150, 200, or 250 *μ*m (*n* = 8). Between the two cells, intercellular collagen density (mg/mL) and AI were measured throughout simulations. Gathering intensity (GI) was calculated by dividing intercellular collagen density by the initial average collagen density (i.e., 1.5, 2.0, or 3.0 mg/mL) [40].

Statistical comparisons were made using graphpad prism 7.0 a (GraphPad Software, Inc., San Diego, CA). Percentage of segments with greater than 1% strain in cell-free simulations and total strain energies in two-cell simulations were analyzed by one-way ANOVA followed by Tukey's correction for multiple comparisons. Results were considered statistically significant for *p* ≤ 0.05.

The model and code are available online by accessing the website.^{3}

## Results and Discussion

### Mechanical Response of Cell-Free Collagen Gels

#### Tensile and Shear Strain.

Our approach to modeling collagen networks is similar to the work by Stein et al. [24] and others [18,25,28]. A collagen gel is treated as a network of individual, cross-linked fibrils. Fibrils were modeled as elastic beams that resist stretching and bending; cross-links were modeled as torsional springs. Thus, comparing results from our model to that of Stein et al. allows us to check our agent-based model's results against well-established finite element methods. Stein et al. observed that simulated gels strain-stiffened consistent with their [24] and others' [47,48] experimental measurements on collagen gels.

In our agent-based model, before strain was applied to a simulated collagen network, all fibrils were perfectly straight and oriented in random directions (Figs. 3(a) and 3(e)). As tensile strain was increasingly applied, the network elongated along the axis of strain and compressed in the perpendicular direction (Figs. 3(b)–3(d)). At the fibril level, two distinct behaviors were observed. (1) Contiguous groups of fibrils that extended across the network successively aligned along the axis of strain. (2) Many other fibrils became compressed and buckled between cross-links. Once fibrils aligned, further tensile strain resulted in stretching of individual fibrils, appearing straight and labeled red when strain exceeds 1% relative to its resting length (*L _{R}*). Though these aligned fibrils appeared to be a single fibril, this could not have been the case because the maximum length for a single fibril was roughly half the width of the undeformed network. Thus, aligned fibrils must be at least two fibrils that are cross-linked. As greater strain was applied to the networks, greater numbers of fibrils appeared to align and become stretched in tension. Networks subjected to shear strain also exhibited fibril alignment and buckling that increased with the level of strain with a few notable differences (Figs. 3(f)–3(h)). Since the

*X*coordinates of nodes on the sheared edges of the network are fixed relative to each other, networks did not collapse as they did when subjected to tensile strain. Therefore, unlike networks subjected to tensile strain, it was possible to see how many fibrils were aligned at high shear strain (Fig. 3(h)). In addition, the initial alignment and stretching of fibrils generally occurred at higher strains than for networks subjected to tensile strain. As collagen concentration was increased from 2.0 to 4.0 mg/mL, the number of fibril segments increased proportionally; however, the percentage of fibril segments experiencing more than 1% strain also increased (Fig. 3(i)).

#### Stress–Strain Relationships.

In response to tensile and shear strain, simulated networks exhibited a biphasic stress–strain curve that is indicative of strain-stiffening (Figs. 4(a) and 4(b)) and qualitatively similar to the response of real and simulated collagen networks reported by Stein et al. (Fig. 4(b)). At low strains, the relationship between stress and strain was linear and elastic modulus increased with collagen concentration (Fig. 5). Stress–strain curves steepened sharply for tensile strains above 0.018 and shear strains above 0.054, similar strains at which fibril alignment, measured by the anisotropy index, began to increase (Figs. 4(c) and 4(d)). New slopes were maintained for strains up to 1.05. For all networks subjected to tensile or shear strain, stress increased with increasing collagen concentration. Removing dead ends on fibrils did not change the equilibrium value of stress for a given strain, but did modestly increase the rate of equilibration (data not shown).

The incremental elastic modulus for a single network subjected to tensile strain is plotted in Fig. 4(e). For strains below 0.05, the modulus was nearly constant, around 0.1 kPa. Between strains of 0.05 and 0.20, the elastic modulus increased almost three orders of magnitude to ∼44 kPa. Above a strain of 0.20, the elastic modulus continued to increase, but more slowly, and reached ∼144 kPa at a strain of 1.05. To help determine the mechanism behind this strain-stiffening behavior, we performed tensile testing of an identical collagen network, but removed the torsional springs (Fig. 4(e)), an idea inspired by Stein et al. [24]. This change meant that nodes along fibrils and at cross-links would act as freely rotating hinges that do not provide resistance to movement. At low strains below 0.05, the network provided negligible resistance to strain and possessed an elastic modulus below 1 *μ*Pa. However, at a strain of 0.10, the elastic modulus increased to nearly match that measured for the same network with torsional springs (3.44 kPa versus 3.77 kPa). For all strains above 0.15, incremental elastic moduli for networks with and without torsional springs were within 1% of each other. Fibril alignment for this network with and without torsional springs was quantified by measuring AI (Fig. 4(f)). AI was constant for both conditions at low strain, but increased progressively for strains above 0.025. This increase in alignment slightly preceded the onset of strain-stiffening. These data indicate that at low strain, stiffness of fibrous networks is predominantly due to fibril resistance to bending and torsion at cross-links. For networks without torsional springs, stiffness can only be attributable to stretching of fibrils. Since a single fibril is modeled to behave as a linearly elastic material in response to tension, strain-stiffening in our model must therefore be due to the progressive alignment of fibrils along the axis of strain and engagement in tension.

In our model, strain-stiffening occurred above a tensile strain of 0.018 and shear strain of 0.054. These are lower than values reported by Stein et al. who found that real and simulated networks began to stiffen above a strain of 0.10 (Fig. 4(b)) [24]. Our networks might exhibit strain-stiffening at lower strains than those studied by Stein et al. since the initial confirmation of our networks consisted of only straight fibrils while the initial conditions of the networks studied by Stein et al. were based on images of actual collagen networks and contained fibrils that were curved. Using finite element models of discrete fiber networks, Onck et al. have shown that decreasing the persistence length of fibers (i.e., increasing fiber curvature) delays the strain at which strain-stiffening occurs, but does not change the overall shape of stress–strain curves [25]. To explore the relationship between fibril curvature and the strain at which strain-stiffening becomes pronounced, we introduced varying degrees of fibril curvature to determine its effect on network mechanical properties. Fibril curvature was introduced to a network with straight fibrils by moving all nodes along fibrils 0.5, 1.0, 1.5, 2.0, and 2.5 *μ*m perpendicular to their fibril's orientation while leaving nodes at cross-links fixed (Fig. 2(d)). Results were averaged from eight networks subjected to tensile and shear strain for each level of curvature. Similar to Onck et al., increasing fibril curvature did not change the shape of stress–strain curves, but did shift the curves to the right, postponing the onset of strain-stiffening. Across the range of fibril curvature tested, we observed a shift in the onset of strain-stiffening from 0.018 to 0.15 in response to tensile strain and 0.054–0.275 in response to shear strain (Figs. 6(a) and 6(b)). These results provide at least a partial explanation for the difference in values of strain at which networks strain-stiffen between our model and Stein's.

Our agent-based model qualitatively predicts the fibril reorganization and strain-stiffening properties of cell-free collagen gels in response to tensile and shear strain. While this is an important validation of our simulated collagen network, it did not assess the ability of our model to handle matrix remodeling due to cell traction. To address this aspect, we compared our model predictions to two additional sets of experimental data.

### Remodeling Around a Single Cell

#### Collagen Displacement.

Mohammadi et al. suspended fluorescent magnetite microbeads (2 *μ*m diameter) within collagen hydrogels cast in nylon frames measuring 1700 *μ*m by 1700 *μ*m [46]. These microbeads acted as fiduciary markers that allowed collagen remodeling to be tracked over the course of their experiment. A single 3T3 fibroblast was then placed onto this collagen network and the displacement of microbeads was tracked at regular intervals. After 6 h, peak bead displacement of 60 *μ*m occurred closest to the cell and displacement decreased with distance up to ∼500 *μ*m from the cell (Fig. 7). The displacement field therefore did not extend to grid boundaries.

We placed a single cell on the center of a free-floating, circular 2.0 mg/mL collagen network 1000 *μ*m in diameter. This network was sufficiently large so that it would exceed the area that could be remodeled by the cell during the simulation. The free-floating network we model is different from the experiments by Mohammadi et al. that utilized constrained collagen networks. However, since the deformation field does not extend to the boundary in either case, we feel that this computational design offers a fair comparison. Instead of tracking beads, we tracked the position of individual nodes on our collagen network. As the cell continuously applied traction force, by repeatedly binding, pulling, and releasing nodes on the network, collagen fibrils steadily accumulated within a region near the cell. After 108,000, 216,000, 324,000, and 432,000 iterations, peak displacement was measured for fibril segments initially at 100 *μ*m (Fig. 7). Due to the large magnitude of fibril displacement close to the cell, on the order of tens of microns, and since most fibril segments close to the cell can only be displaced a maximum of their distance to the cell, displacements of individual fibril segments were not tracked for fibril segments initially within 100 *μ*m of the cell. After 108,000 iterations, collagen was displaced up to ∼225 *μ*m from the cell. The area within which nodes were displaced increased with time beyond a distance of ∼300 *μ*m at 432,000 iterations. Though the magnitude of displacement is greater in the experimental data, the shape of our displacement curve at 432,000 iterations approximates the trend for displacement reported by Mohammadi et al. (Fig. 7). By also plotting data at earlier time points, we were able to observe an additional feature that was not obvious at the 432,000 iterations time point or revealed by the experimental data. Beyond the linear region, our simulations predict displacement will form an asymptotic-like curve and approach zero.

### Remodeling Between Pairs of Cells

#### Gathering Intensity.

McLeod et al. quantified remodeling between pairs of endothelial cells cultured within reconstituted type I collagen hydrogels [40]. Distances between pairs of cells and intercellular gathering intensities were obtained by imaging cell-seeded hydrogels after 4 h of culture. Gathering intensity represents the degree of collagen compaction and is calculated by dividing the average fluorescent intensity in the remodeled intercellular region by the average intensity in a region that has not been remodeled (Fig. 8(a)). By plotting gathering intensities against corresponding cell–cell distances, they observed an inverse relationship between gathering intensity and cell–cell distance (Fig. 8(b)). In addition, they noticed a downward shift in gathering intensity at a constant intercellular distance as collagen concentration was increased in the range between 1.5 and 3.25 mg/mL.

We designed our simulations to match the experimental system used by McLeod et al. Each simulation consisted of two cells within a collagen network (Fig. 8(c)). We varied cell–cell distance and collagen concentration, then quantified gathering intensity after 288,000 iterations. We ran our simulations until we achieved roughly the same magnitude of gathering intensity for a 1.5 mg/mL collagen network. We calculated gathering intensity by quantifying the average concentration of collagen in the intercellular region between cells and dividing this value by the initial average collagen concentration of the network. As with the experimental data, in simulations, gathering intensity appeared inversely related to the initial distance between cells (Fig. 8(d)). We also observed a trend toward a decrease in gathering intensity with increasing collagen concentration. However, for cells separated by less than 100 *μ*m, the differences between measured gathering intensities at different collagen concentrations appeared smaller in magnitude in the computational data than in the experimental data (Figs. 8(b) and 8(d)).

Using the final state of our two cell simulations for 2.0 mg/mL collagen networks, we created heat maps of average strain energy by summing the energy stored in every tensional and torsional spring on each patch (5 *μ*m by 5 *μ*m square), then by averaging the strain energy at each patch for all simulations at each intercellular distance (Fig. 9). By comparing the average strain energy maps for cells separated by different distances, an interesting pattern emerged. For cells separated by 250 *μ*m, the average strain energy appeared roughly uniform around each cell and energy decreased with distance from each cell. As cells were placed closer together, the strain energy field around each cell became increasingly nonuniform and a region of high strain energy in the direction of the other cell increased in prominence, as can be seen by the yellow region overlapping between the two cells (Figs. 9(d) and 9(e)).

The strain energy maps for two nearby cells (e.g., Fig. 9(e)) indicate that within the pericellular region (i.e., <25 *μ*m from the cell surface), areas of lowest strain energy are not on the side of the cell opposite the neighboring cell, as might be expected, but will instead be directly adjacent to the region of peak strain. To better quantify this observation, we graphed the average strain energy 10 *μ*m from the cell membrane nodes as a function of angular position (Figs. 9(f)–9(j)). 10 *μ*m was chosen as the point to sample strain energy because this value approximates the region the simulated cell will interact with when spread and extending pseudopodia. Consistent with the strain energy map, for cells separated by 75 *μ*m, the minimal average strain energy was not in the direction away from the other cell (i.e., 180 deg) but ∼45 deg off of the direction pointing to other cell. The same pattern is also distinguishable for cells separated by 100 *μ*m, but weakens at greater intercellular distances.

We speculated that the minimum strain energy occurred at ∼45 deg due to differences in collagen density as fibrils from this region might be recruited relatively easily to become aligned between the two contractile cells (Fig. 8(c)). To explore this idea, we calculated and mapped the average collagen density on each patch for cells separated by 75 *μ*m and confirmed the density of collagen was relatively low in the areas ∼45 deg from the direction of the other cell compared to other pericellular regions (Fig. S8(*a*), which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection). In contrast to the simulation results, regions of relatively low collagen density adjacent to the aligned collagen are not obvious in the one representative image from the in vitro study (Fig. 8(a)). Since collagen densities in these regions were not quantified, it remains uncertain as to whether this pattern also forms in vitro. The pattern of collagen density observed in simulations suggests that the spatial variations in collagen density after remodeling contribute to the observed distribution of strain energy within the collagen network (Fig. 9). However, the variation in collagen density is only part of the explanation; a map of the average strain energy divided by the average collagen density reveals that the pattern of strain energy stored within individual collagen fibrils also contributes (Fig. S8(*b*), which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection).

### Limitations of the Model.

While our agent-based model reproduces various experimental observations and gives insights into the underlying mechanical basis of these observations, our model also has notable limitations. These limitations can be roughly grouped into three categories: uncertainty in model parameters, failure of the model concepts to accurately represent the real physical system, and computational limitations. Beginning with uncertainty in the model parameters, as noted in the introduction, differences in the conditions during collagen polymerization profoundly affect collagen microarchitecture and mechanics [23,41]. Since precise values for collagen fibril diameter and modulus, two parameters that predominantly determine fibril mechanics, were not known for each experimental condition, we chose to construct networks using a single value for each of these parameters for all simulations, estimated from within experimentally determined ranges. Additionally, the magnitude of cellular traction depends on both cell type and environment [42]. Though traction force microscopy and other methods have quantified total cell traction for various cell types, because of the number of assumptions required, we did not have confidence in translating these data to an appropriate parameter value for the traction force applied by each pseudopodium in our 2D model. Therefore, we empirically chose a value that resulted in a rate of network deformation comparable to experimental data. We could have improved the agreement between our simulation results and a specific experimental data set by judiciously choosing different parameter values. However, since there is no way of knowing whether the specific values we chose more accurately match those from the experimental condition, any improvement in agreement between the simulations and experimental data would not help validate the model.

Another limitation of our model stems from its failure to accurately represent the real physical system, either due to an incomplete understanding of the physical system or intentional simplifications made to facilitate model development. Our incomplete understanding of the physical system extends even to relatively basic concepts. For example, the dynamics of collagen polymerization are not well understood; the experimental literature is unclear as to whether a change in collagen concentration influences fibril diameter [41] and/or fibril number [21]. In addition, fibrils of a given diameter may have different bending moduli depending on the bulk collagen concentration during polymerization [21]. These limitations, founded on an incomplete understanding of the physical system, are shared by other models of cell–matrix interactions. To navigate this incomplete understanding, we chose to construct collagen networks by randomly placing fibrils and model a change in collagen concentration simply as a proportional change in fibril number with a fixed elastic modulus. Many of these assumptions have also been employed by others [18,25,28] and serve as a useful first approximation for simulating collagen networks [18]. Regarding cross-linking, in order for simulations to remain stable, only a portion of fibrils that cross could be cross-linked so as to avoid placing nodes along a fibril too close together. This limitation is discussed in more detail in the Supplemental Material, which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection.

An example of a limitation resulting from an intentional simplification is the “contact problem.” Like others using fiber structural models [18,23–27] and agent-based models of fibers [28], we do not consider the forces that that one fibril would apply on another fibril as they move into contact with one another. As a result, fibril segments can and do cross through one another in our (and others') simulations. There is no fundamental reason why the contact problem cannot be incorporated into the agent-based model, but it would be very computationally intensive, as it is in finite element models. Similarly, we and others have used two dimensions to model an inherently 3D system [18,25,28,31]. While the general agreement between experimental results and our modeling results suggests that using a model with reduced dimensionality can provide insights, a concept reviewed by Mogilner and Odde [49] reduced dimensionality limits in the exploration of some hypotheses. For example, in addition to the contribution of fibril curvature, another possible explanation for the observation that strain-stiffening occurred at a lower strain in our model than in the 3D simulations of Stein et al. is the difference in dimensionality of the networks. Noting that our results and those of others suggest that fibril alignment precedes and contributes to strain-stiffening, one might hypothesize that because fibrils in a 2D network are already in a single plane, less reorientation, and therefore less strain, might be required for fibrils to align along the axis of strain. A model that could simulate both 2D and 3D networks could test this hypothesis, but our model that is currently limited to two dimensions cannot. A discussion of computational consequences of transitioning to a 3D model is provided in the Supplemental Material, which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection.

The uncertainty in model parameters combined with a failure of model concepts to accurately represent the real physical system might explain certain inconsistencies between our computational and corresponding experimental results. Although our model predicted decreasing network compaction between two cells with increasing initial collagen concentration, the magnitude of the effect in the simulation appeared to be less than that seen in the experimental data. This observation may be due to our assumption that a change in collagen concentration results in a proportional change in fibril density. However, we also modeled a change in collagen concentration as a change in fibril diameter, but this approach did not support a greater response by gathering intensity to collagen concentration (Fig. S2, which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection). We show these data and offer other possible explanations for why our model does not match this experimental prediction in the Supplemental Material, which is available under the “Supplemental Data” tab for this paper on the ASME Digital Collection.

### Potential Utility of the Model.

One of the major goals in this work was to validate our agent-based model against various experimental data sets. While validation is an important component for many types of models, validation is especially important for agent-based models since agent-based mechanical models are relatively poorly characterized and few models have been validated against experimental data. Our agent-based model predicted emergent behaviors (e.g., behaviors not intentionally programmed into the model but that arose from the simple mechanical rules specified). Importantly, these emergent behaviors were also evident in the experimental data including (1) strain-stiffening of cell-free collagen networks loaded in response to tensile or shear strain, (2) network stiffness increasing with concentration, (3) fibril displacement toward a single cell with the magnitude of displacement decreasing with distance from the cell, and (4) collagen fibril alignment and compaction in the region between two cells with the degree of compaction decreasing with increasing distance separating the cells. The ability of the agent-based model to predict diverse responses including the organization (e.g., displacement, alignment, and compaction) of fibrils and the mechanical properties of the fibrous networks helps validate the model.

Our model is not unique in its ability to predict some of these emergent behaviors. However, to our knowledge, no other model allows for dynamic and repetitive hand-over-hand traction applied by the cell to the collagen network. By modeling traction as many pulls, our agent-based model cumulatively produced a greater degree of remodeling than is supported by existing finite element models of discrete fiber networks. Accordingly, our agent-based model is well suited for modeling the more biologically realistic, large matrix deformations in the pericellular region observed experimentally (Figs. 7 and 8). For example, in our model of a single cell in a collagen network, collagen fibril displacements exceed 50 *μ*m (Fig. 7) and network remodeling is visually striking (Fig. 8(b)). This value is not a maximum, simply what was achieved during the fixed duration of simulations. In comparison to the relatively large peak displacement achieved by modeling traction as many pulls, one might expect to measure a peak displacement of only ∼10 *μ*m were traction to have been modeled as a single pull.

It should be noted that we used essentially the same model with the same parameters for all simulations with only minor variations to account for specific differences in experimental conditions. Thus, our model was not optimized to fit one particular data set, which suggests that that model predictions might give insights into other analogous physical systems or properties that cannot be measured experimentally, such as strain energy. Mapping the strain energy in the final state of our two-cell simulations shows, as expected, that the strain energy is the greatest in the region between two adjacent cells (Figs. 9(e) and 9(j)). The region with the lowest strain energy is not in the direction away from the other cell (i.e., 180 deg) but ∼45 deg off of the direction pointing to other cell. This pattern of strain energy is interesting in light of the observation that endothelial cells (i.e., the cell type used in the two-cell collagen network remodeling experiment) first realign the fibrils between two adjacent cells before preferentially migrating or elongating toward one another [40,50]. The signals cells detect to direct their migration is not known, but features such as stiffness, binding site density, and alignment are reasonable candidates. Bischofs and Schwarz suggested that cells align so that they exert the maximal work on an elastic substrate [51]. One could consider an analogous case where strain energy, which is directly related to work in a purely elastic solid, where cells migrate in the directions of maximal strain energy density. If cells sense these areas of low strain energy adjacent to the high strain energy region, then their perception of the high strain energy might be enhanced. This pattern might increase the likelihood of migration and persistence along aligned fibrils between the two cells.

## Acknowledgment

The authors thank Hamid Mohammadi, Alisha Sarang-Sieminski, and John Wiley and Sons for permission to include previously published data in this manuscript.

## Funding Data

Division of Civil, Mechanical and Manufacturing Innovation (Grant No. 1334757).