Abstract
Two mathematical models for fibroblast–collagen interaction are proposed which reproduce qualitative features of fibroblastpopulated collagen lattice contraction. Both models are force based and model the cells as individual entities with discrete attachment sites; however, the collagen lattice is modelled differently in each model. In the collagen lattice model, the lattice is more interconnected and formed by triangulating nodes to form the fibrous structure. In the collagen fibre model, the nodes are not triangulated, are less interconnected, and the collagen fibres are modelled as a string of nodes. Both models suggest that the overall increase in stress of the lattice as it contracts is not the cause of the reduced rate of contraction, but that the reduced rate of contraction is due to inactivation of the fibroblasts.
1. Introduction
In 1972, Elsdale and Bard first reported the contraction of collagen gels by fibroblasts [1]. Seven years later, Bell et al. introduced the fibroblastpopulated collagen lattice (FPCL) contraction model [2] with the goal of better understanding how contraction affects closure of an open wound. The FPCL system is intended to mimic cell–matrix interactions that occur in wound granulation tissue. Although much has been learned from this system, there are still many fundamental open questions. In this paper, we develop twodimensional mathematical models designed to investigate the role of cellular forces in contracting lattices.
We describe two mathematical models that predict FPCL contraction for lattices with various cell densities. In §2, we give necessary background including a description of freefloating and constrained lattices, a brief discussion of myofibroblasts and fibroblasts, and a review of some previous modelling efforts. In §3, we present a description of the experimental methods and the numerical software used in the implementation of our model. Section 4 describes two separate mathematical models—one with a collagen lattice structure and one with a collagen string structure. In §5, we describe the results of our numerical simulations, and we conclude with a discussion in §6.
2. Background
Since their introduction, several modifications have been made to FPCLs to model different biologically relevant lattices. The two most frequently studied lattices are freefloating lattices and constrained lattices. Freefloating lattices float on the medium and are allowed to freely deform. Although there are local stresses in these lattices, they are imposed by the collagen structure within the lattice [3]. The constrained lattices have external forces imposed on the lattice which constrain the shape of the lattice. In this manuscript, we restrict our study to freefloating lattices. For a more complete review of FPCLs, the reader is referred to [4].
Fibroblasts in FPCLs exhibit two phenotypes, the normal phenotype (referred to as fibroblast) and the myofibroblast phenotype. The myofibroblast phenotype is characterized by the expression of αSMA, the presence of bundles of contractile microfilaments and extensive celltomatrix attachments. Myofibroblasts appear to exert greater forces than fibroblasts, are more adhesive to the extracellular matrix and therefore less motile, and produce more extracellular matrix proteins [5–7]. The predominant phenotype in a lattice is dependent on experimental design [8].
There are several mathematical models of cell–extracellular matrix interactions that focus on the forces involved. Early models treated the cells and the extracellular matrix as continuum variables and used classical mechanical laws for continuum media to formulate the partial differential equations used in the models [3,9–12]. Later, in an effort to better model the collagen network, Baracos and coworkers [13–15] developed a hybrid method that considers the fibrous structure to determine forces in a volume averaging way and thus deduce the biomechanical properties of collagen tissue. On one scale, they consider the extracellular matrix as a discrete fibrous structure, but in solving the tissue properties, they use a continuum description. Our models, however, do not use a continuum description of the extracellular matrix or the cells.
There are multiple discrete models that are relevant to the model presented in this manuscript. Stein et al. [16,17] model threedimensional collagen structures with discrete fibres. They consider these fibres as stiff rods that can twist at crosslinking points to determine the deformation of the lattice. Schluter et al. [18] takes a more phenomenological approach to model a single cell and discrete fibre interactions to understand how the extracellular matrix affects the migration of cells. In their model, they model a drag force on the cells through the matrix and realign the matrix in the direction the cell is moving. Reinhardt et al. [19] use an agentbased model to simulate both complex extracellular matrix remodelling and durotaxis. They model the force interactions between discrete fibres and cells using the Fruchterman–Reingold algorithm [20] where the links between the cell and the binding site act as springs and the binding sites act as electrically charged particles. This model has the advantage of straightening collagen fibres; however, in this paper, they consider only one or two cells. Finally, there are modelling efforts that use a discrete fibre formulation to derive a closed form for the strain energy [21]. The goals of this last type of work are quite different from ours.
The model presented in this paper is a discrete forcebased model that models the interaction between fibroblasts and a collagen lattice. This model differs from the previously described models as we model an entire lattice with many cells and focus on the contraction of the entire lattice. This model also differs from previously described models as we allow for compaction (a biochemical process) and also for differing cell types (fibroblasts and myofibroblasts).
3. Methods
3.1. Cell lines
Human dermal fibroblast cultures were derived from neonatal foreskin explants and maintained in complete Dulbecco's modification of Eagle's medium (DMEM), with 10% fetal bovine serum (FBS) and 15 μm ml^{−1} of gentamicin. Fibroblasts were studied between their eighth and 12th passage. DMEM and FBS were purchased from Life Technologies (Rockville, MD).
3.2. Castingpopulated collagen lattices
Fibroblast PCL containing either 3000, 10 000, 30 000 or 100 000 cells per ml of 1.25 mg of acid soluble rat tail tendon collagen in 1 mM HCl and complete DMEM. In a 60 mm Petri dish, 0.2 ml drops of cell–collagen–medium mix were pipetted onto the surface of Petri dishes and allowed to polymerize before adding 2 ml of complete DMEM, which freed the lattices from the dish. Between four and five PCLs in each dish were maintained at 37°C, with 5% CO_{2} in a watersaturated atmosphere incubator for 3 days without a change of medium. The diameter of the cell populated collagen lattices was measured initially and every 5 h until 40 h had elapsed.
3.3. Numerical methods
The numerical computations for our models use the following software packages. Initial triangulations for the collagen lattice model are created with Triangle [22]. The system of differential equations used to model the system is solved using SUNDIALS CVODE [23]. To determine the contraction of our simulated collagen, we compare the area of the original lattice with the area of the contracted lattice. These areas are determined by taking the area of the convex hull of the node locations determined by qhull [24]. In order to find model parameters that would give results fitting the data, we used the optimization software gsl [25]. Finally, we used Matlab to visualize our results.
4. Model
Our mathematical model is forcebased, and models both the cells and the collagen lattice. Although there is a significant amount of biochemical remodelling of the extracellular matrix by fibroblasts in both the wound environment and in collagen lattices, in the first 24–48 h after casting an FPCL, there is not much collagen deposition by the fibroblasts. The contraction of the FPCL therefore seems to be due to the forces generated by the fibroblasts.
Although the model is force based, we allow the collagen fibres to compact. The compaction of collagen is the combining of fibrils into longer and thicker collagen fibril. When the distance between two collagen fibrils is small enough, the fibrils will combine into a larger collagen fibre. Eventually, this merging process forms longer and thicker collagen fibres. Evidence suggests this is the fundamental process that underlies wound contraction [26].
4.1. Cells
In our model, we consider a collagen lattice with N cells. A single cell is modelled as K integrinbased adhesion sites (Isites) which exert force on a common location which can be thought of as the centre of mass (figure 1). We note that for the simulations shown in this manuscript we fix the value of K at 10. The Isites exert forces on the centre of mass according to Hooke's law, i.e. the force is proportional to the distance from the centre of mass. We may think of this as if the Isites are attached to the centre of mass with springs which have rest length ℓ, set to be zero except in simulations with the stressdependent contraction mechanism and the stressdependent mechanism (see §5.1). The drag on the centre of mass is modelled as a sphere in water with low Reynolds number and is proportional to its velocity. The equation of motion for the centre of mass of the ith cell is given by 4.1where x_{i} is the location in of the cell centre and represents the velocity of the cell and i = 1, … , N. The spring constant which defines the cell force is α, the rest length of the spring is ℓ and the drag coefficient is C. The righthand side of the equation expresses the forces present on the cell centre from the integrin attachments, and the lefthand side of the equation expresses the motion of the cell. If the Isites remained permanently attached, the cell would eventually reach an equilibrium position and both sides of the equation would equal zero. The Reynolds number is low, thus because of the relative magnitudes of the coefficients, the acceleration term on the righthand side of the equation is set to zero.
We constrain the location of the Isites to be at lattice nodes which are specific locations on the lattice (see §§4.2 and 4.3). The lattice node locations are denoted by y_{k}, and Isites from the same cell or other cells can be attached to the same node location. The set of indices p_{i}_{,1}, p_{i}_{,2}, … , p_{i,K} specify the lattice nodes associated with the Isites of cell i.
The Isites remain attached to a node for a random amount of time, taken from a Poisson distribution with mean 60 s (for simulations with the stressdependent attachment mechanism and the stressdependent mechanism the time is extended, see §5.1). Once an Isite detaches, it immediately reattaches, thus, there are always K attachment sites for each cell. To determine the new attachment node, we determine a direction and a distance from the cell centre. The direction is chosen by choosing a value from a uniform distribution over the interval [−2, 2] and then perturbing the angle of direction for the previous location of the Isite by this value. Once the direction is fixed, the distance is chosen from a uniform distribution between 0 and 115.726 μm, from the cell centre. Once the direction and angle have been determined, we negate the distance with probability 0.2. This has the effect of reversing the direction the Isite reaches in 20% of the cases. Finally, we allow the Isite to attach to the nearest lattice node. For more information about the Isites, the reader is referred to a related model discussed in [27].
4.2. Collagen lattice
The collagen lattice is modelled by several nodes which are connected to form a network of springlike connections (figure 2). To create the collagen lattice, M nodes are placed in the domain, and the connections between the nodes are determined through the use of a Delaunay triangulation [28]. Recall the cell Isites are constrained to be at lattice nodes, as shown in figure 3.
Forces acting on lattice nodes come from the lattice or the cells. The equation of motion for the lattice node k is 4.2The first summation on the righthand side denotes the forces owing to connections to other lattice nodes. Only the nodes directly connected to node k can have nonzero forces. The second summation on the righthand side gives the forces owing to cells which have Isites bound to node k (figure 3).
Lattice forces come from connections with other nodes and are classified into two types: normal or compacted. Having the second type of connection allows the compaction of collagen and is a nonreversible process. The forces owing to both types of connections are springlike in that if the connection is stretched the force is proportional to the stretching. If a normal connection is compressed however, there is no force generated. This assumption is due to the nature of collagen. When the collagen fibrils are pulled, they resist the pulling owing to their association with other fibrils. Yet, if a cell exerts forces at two points along the same fibril drawing the two points closer, the fibril is not compressed, but becomes slack between the two points similar to a rope. Compacted links act as true springs, in that, when a link is compressed, a force is exerted. Thus, the force owing to a lattice connection between node k and node m is defined as
4.3
Here, y_{m} is the location of lattice node m, is the distance between node k and node m, ℓ_{k,m} is the rest length of the connection between node k and node m and is set as the initial distance between the nodes at the beginning of the simulation, β is the spring constant for normal links, β^{*} = d_{β}β is the spring constant for compacted links, and is the rest length for the spring connecting node k with node m when the link is compacted. Initially, all links are normal and become compacted if the distance between two linked nodes becomes small enough, that is, if . When links are compacted, the rest length of the spring is shortened (d_{ℓ} < 1), the spring constant is increased (d_{β} > 1) and the link resists compression.
The forces due to the cell i are defined by 4.4where x_{i} is the cell centre location, y_{pi,j} is a lattice node location, α is the spring constant and ℓ is the rest length of the integrin. Here, δ(0) = 1 and δ(x) = 0 for any nonzero x and indicates whether the jth Isite of cell i is interacting with node k.
4.3. Collagen strings
We now introduce an alternate model formulation for the collagen lattice. Rather than place M nodes arbitrarily in a domain in , and then interconnect them using the Triangle program we instead create collagen strings. These strings can be thought of as a chain of collagen nodes which are connected by normal connections, that is, the connections behave as springs when stretched but do not resist compression (see description of normal connections in §4.2).
To create a collagen string lattice, we first fix a distance between the nodes of 100 μm and set a string length by specifying the number of lattice nodes in the string. We have the option of arranging the strings in an orderly fashion to mimic spun collagen [29]; however, to better compare this model to the lattice model, we instead allow them to arrange in a random fashion. We do this by randomly selecting a starting position in our domain and selecting the angle formed between two segments of our string from a uniform distribution in the interval [π/6, 11π/6] to avoid kinks. We continue this process until a string of the desired length is created. We repeat the process for each collagen string until the desired density of collagen nodes is reached.
We now note a fundamental difference between the collagen string model and the collagen lattice model. In the lattice model, all node linking is done before any simulations are run. In the collagen string model, nodes are allowed to link at each time step if they are within 5 μm of each other. When nodes are linked in this scenario, they form a compacted connection.
4.4. Model parameters
Our model contains many parameters, but the most important parameters for accurately modelling collagen lattice contractions are

— α which represents the force the cells exert on the lattice,

— β which is Young's modulus for normal collagen, and

— β^{*} which is Young's modulus for compacted collagen.
These parameters have a large range of values reported in the literature.
The force that fibroblasts exert on the extracellular matrix (i.e. the value of α) is measured using four different techniques [7]. Each technique gives very different values which results in a wide range of reported fibroblast forces in the literature from 1 nN to 2.65 μN. The first technique is to measure the deformation of a silicon substrate by fibroblasts [30,31]. This technique gives the highest force measurements. The second method is to measure the deformation of micromachined devices by fibroblasts [32,33]. The third technique measures forces on FPCLs and determines the average force exerted by a single fibroblast [34,35]. The previous two techniques report forces in the range of 0.1–138 nN. The fourth technique uses column buckling theory to determine the force of individual fibroblasts in fabricated lattices [36]. They find that fibroblasts exert average forces ranging from 11 to 41 nN with an upper limit of 450 nN. With such a wide range of values reported in the literature, we choose to optimize the parameter α with others to try to best fit the experimental data.
Young's modulus for collagen gels (which depend on the makeup of the gel and whether the modulus is compressive or tensile) ranges from 0.004 to 24 kPa [3,37–41]. Like the spring constant for the cells, this parameter will be optimized in order to fit experimental data.
We chose a value of 60 s for the average duration of the Isite attachments arbitrarily. The only data for this parameter which we are aware of are for Dictyostelium cells and cervical cancer cells [42,43]. We ran simulations with values varying from 6 min to 10 s for the average duration of the Isites. Although the contraction graphs varied, the qualitative conclusions remain valid for all values used.
5. Results
The first goal of our work is the determination of appropriate numerical parameters for α, β and β^{*}, so that the numerical simulations will match the experimental observations. Our experimental method is detailed in §3, and the results for lattices with 3750, 10 000, 30 000 and 100 000 cells per ml, gathered over a period of 40 h, are shown in figure 4. For the numerical simulations, we assumed that only fibroblasts exist in the collagen lattice. We ran our numerical simulations for contraction with both collagen lattices and collagen strings with four different cell densities matching those in the experimental data.
5.1. Collagen lattice results
Scaling the equations by 1/β and using leastsquares to determine the best fit we varied the following values to find the best fit to the experimental observations: C/β, C is the cell drag coefficient; γ/β, γ is the collagen drag coefficient; d_{ℓ} the factor by which the collagen rest length is decreased when it is compacted; d_{p} the factor which determines when the links compact; α/β, α is the spring constant determining the force the cell exerts on the collagen and d_{β} the factor which increases the spring constant of the collagen when it is compacted. Although we were able to find parameter values that gave a good fit to any single cell density, we were unable to find parameter values that were valid for more than one cell density. We refer to these simulations as standard simulations.
This inability to find a set of parameters valid for all the cell densities motivated a closer study of the experimental results. In these results, there are two main features: an initial fast decrease in size of the lattices and then a levelling off as the lattice size stabilizes and decreases less quickly. The main difficulty with our initial optimization was matching both these behaviours. In order to match both features, we altered the model to allow the cells to have stressdependent Isites and contraction.
We modified our model to simulate three other cases. The first modification allowed the attach time of an Isite to depend on the force applied to the Isite and is referred to as the stressdepended attachment mechanism. Although this is not experimentally verified, owing to the mechanosensing ability of cells, it is a reasonable assumption. For the sake of simplicity, we set the attach time to be longer than the simulations effectively making the Isite permanently attached. To implement this modification, we introduce two new parameters to the model. The first, A_{l}, is the minimum length the integrin must be stretched, and the second A_{f} is a force factor. Both components are necessary, because the forces on a lattice node depend on the lattice forces owing to entanglement and the force owing to the cells. For example, if the Isite is attached to a highly entangled region of the lattice and the Isite is far from the cell centre, the net force on the lattice node may be low, yet if the Isite is close to the cell centre and the lattice not entangled, the net force would still be low. Thus, a low net force for a distant Isite indicates a region of high entanglement in the lattice, but a low net force for a close Isite may not indicate a region of high lattice entanglement. The precise rule for the stressdependent attachment mechanism is: if the net force on the node (as determined by its velocity) is less than A_{f}A_{l}α and where i denotes the cell and j denotes the Isite, then the Isite becomes permanently attached. Implementing this rule and optimizing on the parameters A_{f} and A_{l} in addition to the original parameters still did not give the desired results. Cells with the stressdependent attachment mechanism did not simulate the data any better than the standard simulations.
In the next set of simulations, we implemented a stressdependent cell contraction mechanism. This modification was made by allowing the cells to become stronger when the same threshold as the stressdependent attachment mechanism to fix the Isites was reached. In simulations with the stressdependent cell contraction mechanism, the Isites were not fixed, but the cells' contraction properties changed so that they exerted greater forces and resisted compression. This was done by multiplying the force the cells exert, α, by 10 and by changing the rest length of the Isite adhesion to be its length at the first time the rule indicated stronger forces. The simulations with cells having the stressdependent cell contraction mechanism did not give satisfactory results. Simulations that did not fix the length of the lsites to be the new length but had the greater force also failed.
Our final simulation combined the stressdependent attachment mechanism and the stressdependent cell contraction mechanism to create the stressdependent mechanism. Optimizing the parameters where cells had both the stressdependent attachment mechanism and the stressdependent contraction mechanism gave results that matched the data (figure 4 solid lines). Figure 4 also shows the effect of removing either or both of the new cell mechanisms. The simulations in figure 4 use the same parameter set which was optimized for the stressdependent mechanism simulations. These results suggest that allowing cells to resist compression and have permanently attached Isites is important in matching the rate of cell contraction in the freefloating lattices.
Having found simulations that matched the data, we examined the cell distribution at the conclusion of the simulations using the stressdependent mechanism. Figure 5 shows the initial cell distribution in the left panels for each cell density and the final cell distribution in the right panels for each cell density. Figure 6 shows the results from a simulation with 500 000 cells per ml. After 40 h, the cells are much more concentrated near the boundary of the lattice with the cell density in the interior of the lattice being about 84 000 cells per cm^{2}, and the cell density at the edge being about 250 000 cells per cm^{2}. The lattice is also compacted in a tight ring around the edge. Both these features are seen experimentally [3,44]. The experiments by Ehrlich et al. use a cell density of 500 000 cells per ml, and the experiments of Simons et al. have a cell density of 50 000 cells per ml. To determine whether the same trend held for the simulations with 100 000 cells per ml, we averaged the interior cell density and the boundary cell density for 102 runs with the same initial lattice configuration, but different random initial positions for the cells and different instantiations for the random behaviour of the cell motion. The results confirm the same characteristic of the final cell distribution. The average cell density at the boundary of the lattice at 40 h was 31 423 cells per cm^{2}, and the average cell density for the interior was 20 074 cells per cm^{2}. Figure 7 is a magnification of the edge of the lattice shown in figure 5 bottom right panel. Again, a dense ring of collagen can be seen at the edge, just not as prominent as in the higher density case.
We also examined the amount of force each cell was exerting at the final time of the simulation (40 h) and the spatial distribution of the cells. The results are shown in figure 8. The insets show the scatter plot of the cell forces plotted against the radius of the lattice. The average force the cells exert with respect to radius seems to remain relatively constant. This should not be confused with the stiffness of the lattice which likely varies. The contraction data indicate that the lattice is approaching an equilibrium of forces (i.e. the cell forces and the lattice forces are equilibrating and the contraction is levelling off).
5.2. Collagen string results
Once a set of valid parameters were found for the collagen lattice, we wished to simulate the same four cell densities using the collagen string model. We began by simulating the standard case with the expectation that again a poor fit would be found to the experimental data. As can be seen in figure 9, the standard case results in excessive contraction. It is worth noting though that when the collagen string case is compared with the collagen lattice case less contraction occurs for the string case. The reason for this decrease in contraction is due to the fact that the collagen strings are less interconnected initially. As the strings are contracted compaction occurs and the strings become interconnected, however, there still exist collagen nodes that are connected to only one or two other nodes. This creates ‘collagen fingers’ which trail the contraction occurring in the model. For a visualization of these fingers, see figure 10.
Knowing that the stressdependent mechanism was critical to good experimental match for the collagen lattice model we ran simulations using the stressdependent mechanism. Our goal was to minimize the number of parameters that must be changed in order to match the experimental data. We discovered that it was possible to closely match the experimental data by changing only the force factor A_{f}. Recall that the force factor was used in the stressdependent rule to change the cell behaviour, that is, if the net force on the node (as determined by its velocity) is less than A_{f}A_{l}α and where i denotes the cell and j denotes the Isite, then the Isite becomes permanently attached. The values of A_{f} used were 0.00286718, 0.005818, 0.026718 and 0.126718 for cell densities of 3750, 10 000, 30 000 and 100 000 cells per ml, respectively. In order to maintain the closest possible match to the parameters used for the collagen lattice case, we allowed for different values of A_{f}. The reason for the higher force factor for the higher concentration is that with so many cells pulling on collagen nodes a higher tension exists initially. The higher value of A_{f} allows for cells to contract the collagen lattice for a short time until the tension surpasses a factor of the original tension. In the lower density cases, the initial tension is lower, so a lower force must be used to trigger the conversion to permanent Isites and cells which resist compression and are stronger. It is worth noting that for the collagen lattice the force factor A_{f} is 0.172 which is greater than any of the force factors used for the collagen strings indicating a higher initial tension. This is to be expected, because the collagen lattice is more highly connected and moving one node will typically result in moving many nodes requiring a greater force. The greater connectivity of the collagen lattice distributes the cell and collagen forces over a larger region of the lattice.
We also wished to simulate cells with the stressdependent attachment mechanism and the stressdependent contraction mechanism. For simulations with the stressdependent attachment mechanism, we decreased the motility of the cells once our stress rule was met by setting the attach time longer then the simulation. This case is shown by the dashed lines in figure 9. As shown in figure 9 for the collagen strings, the fixing of the integrins seems to be the more important property. To run simulations with cells which have the stressdependent contraction mechanism, we strengthen the cell and allowed them to resist compression (but did not permanently attach the integrin) when the stress rule was met. The results of this simulation can be seen in figure 9. Although not as important as the fixing of the integrins, strengthening the integrins and having the cell resist compression results in significant improvements in experimental matching, especially as the density of the cells increases. Figure 11 shows the initial and final configuration of the lattices for simulations with varying cell density.
6. Discussion
The mathematical models presented here offer new insights into how fibroblasts contract collagen lattices that are consistent with previous modelling work. Before performing the simulations, it was unclear whether the normal tension within a lattice as it contracted was sufficient to explain the experimentally observed rate of contraction. Using two different formulations for the collagen lattice, this work shows that the cells have sufficient force to contract the lattices to a greater degree than what is seen experimentally. The cells must for some reason become less active during the course of lattice contraction. While we cannot prove the necessity of the transition from active cells to inactive cells, the collagen lattice model indicates that it is difficult to explain the behaviour of the collagen without allowing the cells to become both immobile and stronger in resisting both stretching forces and compressional forces. Experiments show that the first 5–10 h for all but the lowest density lattices have the greatest rate of contraction. After this initial period, the lattices contract at a slower rate. The signal causing the slowing rate of contraction may be the greater local tension in the lattice. This may signal the cells to become less active in contracting the lattice, for example by becoming myofibroblasts or some other less actively contracting phenotype. These conclusions are consistent with the results of Stevenson et al. [45], who proposed that a cell in the lattice contracts only the collagen local to the cell and once nearby collagen is sufficiently compacted the cell becomes dormant. Using these assumptions, they were able to accurately model the cellmediated contraction of lattices (not time course data) with different collagen and cell densities. It is also consistent with the modelling results of Reinhardt et al. [19] which suggest that the cells compact the matrix in the pericellular region to a greater extent than other areas. We suggest that as fibroblasts contract the lattice the local collagen structure becomes denser and less easily deformed. The cells in essence are in a stiffer, less compliant part of the matrix (which they formed) and begin to become less mobile. The new behaviour of the cells results in a slowing of the rate of contraction of the entire lattice. Because the myofibroblasts phenotype is more adhesive, less motile and stronger than the normal fibroblast phenotype, one way to test this theory would be to collect timedependent data for the contraction of freefloating lattices where the fibroblasts–myofibroblast transition is promoted (perhaps by adding TGFβ). We hypothesize that a freefloating lattice populated predominantly by myofibroblasts would have a slower rate of contraction.
It is interesting to note the different results in the two model formulations. The collagen lattice model gives an initial lattice that is more highly interconnected than the initial state of the collagen in the collagen string model. In the two higher density cases, the collagen string model shows that the rate of contraction slows down dramatically before 5 h. It is this change in rate which allows the good fit to the data. By either fixing the Isites or giving them a nonzero rest length, the contraction of the lattice slows down dramatically. The same dramatic changes are not seen in the collagen lattice simulations. Of course, the forces are distributed over a larger region in the collagen which should smooth out the transition.
The collagen lattice model not only matches the time course of freefloating FPCLs for several densities, but also show density features observed in experimental lattices. The mathematical simulations with higher cell densities show that the final cell density is higher near the lattice boundary and the collagen is highly organized in a ring surrounding the lattice. It appears that the free boundary is more readily contracted. An explanation for this is that the cells in the centre are anchored by the surrounding cells and the surrounding collagen. In the centre, cells are pulling on collagen nodes in each direction, but these collagen nodes are being pulled by other cells or other collagen nodes resulting in very little movement. Cells on the edges also pull on collagen nodes however those nodes on the edges are ‘free’ from other cells and other collagen connections and contract into the cell. This results in a net movement of these cells towards the centre. After a sufficient number of time steps, this results in the increased density that can be seen.
This theory can be supported by the increased density seen in the centre of the collagen string. Here, the collagen nodes are more free because they are less interconnected. Further evidence supporting this idea is the more highly stretched collagen in the very high cell density case right after the collagen ring on the edge (figure 6b) when compared with the same region in the lower cell density case (figure 7). There are enough cells at the periphery of the high density case to cause the greater stretching of the interior lattice. If the free edge is the cause, an annulusshaped lattice should show the same cell density and collagen rings on both the outer edge and the inner edge.
In summary, both the collagen lattice and the collagen spring model show promise in mathematically modelling the contraction of collagen. An important characteristic of these models is the ability to change behaviour of the cells to become stressdependent. These results show the potential for simple mathematical models to provide insight into biological processes where cells interact with an extracellular matrix.
 Received June 4, 2014.
 Accepted July 24, 2014.
 © 2014 The Author(s) Published by the Royal Society. All rights reserved.