A model for one-dimensional morphoelasticity and its application to fibroblast-populated collagen lattices

TL;DR

This paper introduces a novel one-dimensional morphoelastic model for biological tissue remodelling, capturing structural changes and dynamics such as contraction and re-expansion in collagen lattices, bridging classical mechanics and plasticity theory.

Contribution

It presents a new morphoelastic framework based on deformation gradient decomposition, modeling tissue remodelling from discrete cells to continuum, applicable to collagen lattice behaviour.

Findings

Successfully models collagen lattice contraction and re-expansion.

Bridges classical elasticity with plasticity theory for tissue remodelling.

Provides qualitative agreement with observed biological tissue dynamics.

Abstract

The mechanical behaviour of solid biological tissues has long been described using models based on classical continuum mechanics. However, the classical continuum theories of elasticity and viscoelasticity cannot easily capture the continual remodelling and associated structural changes of biological tissues. Furthermore, models drawn from plasticity theory are difficult to apply and interpret in this context, where there is no equivalent of a yield stress or flow rule. In this work, we describe a novel one-dimensional mathematical model of tissue remodelling based on the multiplicative decomposition of the deformation gradient. We express the mechanical effects of remodelling as an evolution equation for the 'effective strain', a measure of the difference between the current state and a hypothetical mechanically-relaxed state of the tissue. This morphoelastic model combines the simplicity and interpretability of classical viscoelastic models with the versatility of plasticity theory. A novel feature of our model is that while most models describe growth as a continuous quantity, here we begin with discrete cells and develop a continuum representation of lattice remodelling based on an appropriate limit of the behaviour of discrete cells. To demonstrate the utility of our approach, we use this framework to capture qualitative aspects of the continual remodelling observed in fibroblast-populated collagen lattices, in particular its contraction and its subsequent sudden re-expansion when remodelling is interrupted.

Tables

Table 1. Table 1: Summary of the experimental values for the number of fibroblast cells in the collagen gel and the initial gel diameter, as well as the estimated values of the initial cell traction stress σ 0 subscript 𝜎 0 \sigma_{0} and the parameters μ ¯ ¯ 𝜇 \bar{\mu} , τ ¯ ¯ 𝜏 \bar{\tau} and e ¯ crit subscript ¯ 𝑒 crit \bar{e}_{\rm crit} , assuming ℰ 0 = 1000 subscript ℰ 0 1000 \mathcal{E}_{0}=1000 Pa and F cell = 500 ​ n subscript 𝐹 cell 500 𝑛 F_{\text{cell}}=500n N/cell. The remaining model parameters θ ¯ ¯ 𝜃 \bar{\theta} and k ¯ ¯ 𝑘 \bar{k} are optimized to fit the model to each set of experimental data (see Table  2 ).

Experiment	C (cells)	d (cm)	${\sigma }_0$ (Pa)	$\overline{\mu }$	$\overline{\tau }$	${\overline{e}}_{\mathrm{crit}}(\times {10}^{-3})$
Bell et al. (1979)	$7.5\times {10}^6$	$5.3$	$1700$	$2.78$	$1.699$	$16.99$
Talas et al. (1997)	$2.5\times {10}^5$	$2.18$	$334.9$	$2.78$	$0.3349$	$3.35$
Feng et al. (2003)	$2.5\times {10}^6$	$10$	$159.2$	$2.78$	$0.1592$	$1.59$
Guidry & Grinnell (1985)	$1.0\times {10}^5$	$1.2$	$442.1$	$2.78$	$0.4421$	$4.42$

Table 2. Table 2: Values of the model parameters θ ¯ ¯ 𝜃 \bar{\theta} and k ¯ ¯ 𝑘 \bar{k} , obtained from the best fit of the results of the simulations to four different sets of experimental data. Goodness of fit tests are performed by comparing simulations performed using these parameter values with the corresponding experimental data and the resulting χ 2 superscript 𝜒 2 \chi^{2} values for each case are displayed. In addition, we display χ 2 superscript 𝜒 2 \chi^{2} values for the cases where one of θ ¯ ¯ 𝜃 \bar{\theta} and k ¯ ¯ 𝑘 \bar{k} are scaled by factors of ± 10 % plus-or-minus percent 10 \pm 10\% (indicated by a corresponding superscript ± plus-or-minus \pm ), to characterize the robustness of the individual fits.

Experiment	$\overline{\theta }$	$\overline{k}$	${\chi }^2(\overline{\theta },\overline{k})$	${\chi }^2(\overline{\theta },{\overline{k}}^{-})$	${\chi }^2(\overline{\theta },{\overline{k}}^{+})$	${\chi }^2({\overline{\theta }}^{-},\overline{k})$	${\chi }^2(\overline{{\theta }^{+}},\overline{k})$
Bell et al. (1979)	$0.89$	$106.62$	0.0147	0.0313	0.0298	0.0244	0.0218
Talas et al. (1997)	$5.10$	$1294.92$	0.0037	0.0071	0.0063	0.0045	0.0044
Feng et al. (2003)	$16.14$	$183.85$	0.0337	0.0446	0.0434	0.0363	0.0356
Guidry & Grinnell (1985)	$0.33$	$2.97$	0.0590	0.0704	0.0685	0.0921	0.0890

Equations

$$\mathbf{F}=\mathbf{A}\mathbf{G}.$$

$$\frac{\partial e^E}{\partial t}+\frac{\partial }{\partial x}\left(e^{E}v\right)=\frac{\partial v}{\partial x}-g(x,t),$$

$$\sigma =\mathcal{E}(\rho )e^E+\mu \frac{\partial v}{\partial x},$$

$$e_{\text{avg}}=\frac{l_u-l_z}{l_u}.$$

$$\frac{\partial {\sigma }_u}{\partial x}={\tau }_c\left(\delta (x-l_c)-\delta (x+l_c)\right),$$

$${\sigma }_u=\{\begin{array}{ll}{\sigma }_{\text{out}},& l_c<|x|<l_u,\\ {\sigma }_{\text{out}}-{\tau }_c,& |x|<l_c.\end{array}$$

$${\sigma }_{\text{out}}=\mathcal{E}(\rho )e_{\text{out}}+\mu \frac{\partial v}{\partial x}.$$

$${\tau }_c={\sigma }_{\text{out}}-{\sigma }_{\text{cells}},$$

$$e_{\text{cells}}=\frac{{\sigma }_{\text{out}}-{\tau }_c}{\mathcal{E}(\rho )}.$$

$$g(x,t)=-\theta n{\left(-\frac{{\sigma }_{\text{out}}-{\tau }_c}{\mathcal{E}(\rho )}-{\hat{e}}_{\text{crit}}\right)}^{+},$$

$${\hat{e}}_{\text{crit}}<-\frac{{\sigma }_{\text{out}}-{\tau }_c}{\mathcal{E}(\rho )}$$

$$l_u-l_z=e_{\text{cells}}{l}_c+e_{\text{out}}(l_u-l_c).$$

$$e_{\text{avg}}=e_{\text{out}}+\frac{l_{c}}{l_{u}}\big{(}e_{\text{cells}}-e_{\text{out}}\big{)}.$$

$$e_{\text{out}}=e_{\text{avg}}+\frac{n{\sigma }_c}{\mathcal{E}(\rho )}.$$

$${\sigma }_{\text{out}}=\mathcal{E}(\rho )e_{\text{avg}}+\mu \frac{\partial v}{\partial x}+n{\sigma }_c.$$

$$\frac{\partial e}{\partial t}+\frac{\partial }{\partial x}\left(ev\right)=\frac{\partial v}{\partial x}+\theta n\left(-\frac{\sigma -n_l{\sigma }_c}{\mathcal{E}(\rho )}-{\hat{e}}_{\text{crit}}\right),$$

$$\sigma =\mathcal{E}(\rho )e+\mu \frac{\partial v}{\partial x}+n{\sigma }_c,$$

$$\frac{\partial \sigma }{\partial x}=0,$$

$$v(x,t)=\frac{\mathrm{D}\mathrm{u}}{\mathrm{D}\mathrm{t}}=\frac{\partial u}{\partial t}+v\frac{\partial u}{\partial x},$$

$$n(x,t)$$

$$\rho (x,t)$$

$$\mathcal{E}(\rho )={\mathcal{E}}_0+k\frac{\rho -{\rho }_0}{{\rho }_0},$$

$$\sigma (l(t),t)=0,$$

$${\int }_0^{l(t)}w(\xi ,t)\mathrm{d}\xi =u(l(t),t)=l(t)-l_0.$$

$$l(t)=\frac{l_0}{1-w}.$$

$$\frac{\partial w}{\partial t}+\frac{\partial }{\partial x}\left(wv\right)=\frac{\partial v}{\partial x}.$$

$$x^{\ast }=\frac{x}{l_0},$$

$$e^{\ast }=e,$$

$$\bar{\mu }=\frac{\mu }{{\mathcal{E}}_{0}h},$$

$$\bar{\theta }=h\theta n_l,$$

$$U(X,T)=x(X,T)-X,V(X,T)=\frac{\partial U}{\partial T}.$$

$$\frac{\partial }{\partial x}$$

$$\frac{\partial }{\partial t}$$

Full text

††journal: Biomechanics and Modeling in Mechanobiology

1 Introduction

In this study, we present a one-dimensional mathematical model of biological tissue remodelling, based on the multiplicative decomposition of the deformation gradient. An important feature of our model is that the mechanical effects of remodelling are expressed in terms of an evolution equation for the ‘effective strain’ – a measure of the difference between the current state and a hypothetical mechanically-relaxed state of the tissue. This morphoelastic model combines the simplicity and interpretability of classical viscoelastic models with the versatility of plasticity theory. To demonstrate its utility, we show that this model can quantitatively capture aspects of the mechanical behaviour of fibroblast-populated collagen lattices reported in previous experiments.

As biological tissues deform continuously when subjected to mechanical forces, their physical behaviour is often modelled using classical continuum mechanics (Murray, 2001). However, unlike classical solids, living tissues may contain cells that can modify the fundamental mechanical properties of their physical environment. There are a number of processes, most notably tissue growth, in which cells cooperatively alter the tissue structure, changing the relationship between stress and deformation (Chen and Hoger, 2000). Indeed, many tissues undergo a continual process of internal revision and mechanical restructuring, often referred to as ‘remodelling’ (Taber, 1995), in which physical properties of the material, including anisotropy and stiffness, evolve over time. This active remodelling is thought to be significant in a wide range of biological processes such as embryo development and morphogenesis (see Patwari and Lee (2008) for a recent discussion). In particular, it is well known that fibroblast cells, which are found in the stroma of numerous tissues, actively remodel the surrounding extracellular matrix (ECM) by synthesising and reorganising collagen fibres, and that this remodelling is essential for tissue homeostasis and for wound repair (Grinnell, 2003; Majno and Joris, 2004). Remodelling also affects the mechanical stresses experienced by cells in a tissue, which can subsequently modify aspects of cell behaviour. For instance, fibroblasts are known to change their morphology (Gabbiani, 2003; Tamariz and Grinnell, 2002; Tomasek et al, 2002) and phenotype (Amadeu et al, 2003; Gabbiani, 2003; Tomasek et al, 2002) in response to external mechanical cues. Most significantly, the mechanical stresses experienced by fibroblasts during the wound healing process stimulate them to differentiate into more contractile forms: protomyofibroblasts and myofibroblasts (Desmoulière et al, 2005; Gabbiani et al, 1972; Tomasek et al, 2002), which play a major rôle in the subsequent contraction of the wound. The interplay between mechanical stress and active remodelling is hence critically important in wound healing, as excessive contraction is known to lead to pathologies, such as hypertrophic scars and contractures (Roseborough et al, 2004; Enoch and Leaper, 2005; Murphy et al, 2011b).

Since the first descriptions of the kinematics of biological growth the 70s and early 80s (Taber, 1995; Humphrey, 2003), several theoretical frameworks have been proposed to capture the dynamics of remodelling. A very significant advance in this direction was made in the mid 1990s, when Rodriguez et al (1994) and Cook (1995) independently developed a mathematical framework for remodelling that utilizes the multiplicative decomposition of the deformation gradient – an idea that dates back to the 50’s (Bilby et al, 1957; Kröner, 1958, 1959), and which was formalized in the 60’s by Stojanović et al (1964, 1970) and Lee (1969). In this paper, we use the notation developed by Goriely and coworkers (Goriely and Ben Amar, 2007; Goriely et al, 2008; Vandiver and Goriely, 2009; Goriely and Moulton, 2011), where $\boldsymbol{A}$ and $\boldsymbol{G}$ represent the elastic and plastic/growth parts of the deformation gradient, respectively. The decomposition of the deformation gradient tensor thus yields the expression

[TABLE]

The physical interpretation of this decomposition is depicted in Figure 1: $\boldsymbol{G}$ is the deformation gradient tensor associated with a hypothetical deformation from the fixed reference state to a state where all internal stresses are relieved, while $\boldsymbol{A}$ represents the deformation from this ‘zero stress state’ to the current state. For a purely elastic material, the fixed reference state will also be the zero stress state, and we find that $\boldsymbol{G}\equiv\boldsymbol{I}$. However, plastic flow or remodelling enables the zero stress state, and hence $\boldsymbol{G}$, to evolve. A detailed discussion of the fundamental concepts that underlie the multiplicative decomposition of the deformation gradient, is presented in Appendix A.

This approach not only provides a clear and coherent way of understanding growth, but also leads to a natural way of describing “residual stress” – stress that persists even when all loads are removed. Although Fung (1993) had noted much earlier that some tissues (such as arteries) are structured so that it is impossible for the entire tissue to be free of residual stresses unless cuts are made, limited attempts had been made to describe these stresses mathematically. An alternative paradigm was suggested by Goriely and Ben Amar (2007), who coined the term morphoelasticity to describe the combination of elastic and “plastic” changes that are the result of biological growth and remodelling. Despite some similarities, morphoelasticity is quite distinct from classical plasticity. For example, morphoelastic remodelling will generally occur throughout a tissue, not just in those regions where a yield stress is exceeded. Moreover, morphoelasticity can involve changes to the total mass of a tissue (in tissue growth, for example) and/or increases in internal energy, while plastic flow is always mass-conserving and dissipative. This framework has been recently used to mathematically describe aspects of wound healing, namely wound contraction and scar formation (Yang et al, 2013) and dermal wound closure (Bowden et al, 2016), as well as growth in other biological tissues such as axons (García-Grajales et al, 2016).

While remodelling is important in a wide range of phenomena in vivo, surprisingly few in vitro experiments have been developed to explicitly study the macroscopic consequences of remodelling. A notable counter-example is the investigation of the contraction of fibroblast-populated collagen lattices (FPCLs) – cultured fibroblasts embedded in (or placed on top of) three-dimensional ($3$-D) collagen matrices. It has long been known that FPCLs can contract to a small fraction of their initial size within a few days (Bell et al, 1979). This contraction was first observed to be ‘permanent’ by Grinnell and coworkers (Grinnell and Lamke, 1984; Guidry and Grinnell, 1985, 1986), who validated this finding by performing several experiments with varying numbers of fibroblast cells that were placed on top of collagen lattices of different initial densities. Most fibroblasts did not invade the lattice, and were instead found to spread over its surface while reorganising proximal collagen fibres in the direction of spreading. It was thus proposed that the reorganisation of the lattice away from the cells was chemically mediated by the secretion of cell-binding factors such as fibronectin and proteoglycans (Grinnell and Lamke, 1984). It was also observed that the addition of cytochalasin D, which suppresses gel reorganisation by inhibiting cell motility, resulted in a partial re-expansion of the gel (Guidry and Grinnell, 1985, 1986). The relative magnitude of the re-expansion was found to be smaller in gels that were contracted by fibroblasts for a greater period of time, which suggested that the collagen gels were first physically reorganised by fibroblasts and then stabilised by the continued presence of these cells.

In this work, we develop a mathematical description of the contraction of FPCLs that incorporates the mechanical effects of remodelling. In Sec. 2 we provide a detailed summary of FPCLs and the different mechanisms by which they can permanently contract. Then, in Sec. 3, we construct an expression for the rate of morphoelastic contraction of an FPCL, based on a plausible microscopic mechanism of cells rearranging the fibres of the lattice. On varying the two free parameters of this simplified model, we quantitatively replicate features of previous experiments on contracting FPCLs. Finally, in Sec. 4 we discuss how our framework is significantly advantageous in comparison to many other approaches to describing remodelling in biological tissues, and detail possible extensions to our model.

2 Fibroblast-populated collagen lattices

2.1 History and classification

FPCLs were developed by Elsdale and Bard (1972), who used them as a means of investigating fibroblast behaviour in a setting that closely resembles their natural environment. A scanning electron microscope image of a fibroblast embedded in a collagen lattice, taken from Rhee and Grinnell (2007), is shown in Fig. 2. Experiments on such lattices can provide insight into the mechanical interactions between fibroblasts and surrounding collagen fibrils. Since first being developed, these lattices have been used to study the traction forces exerted by fibroblasts in mechanically-relaxed (Bell et al, 1979; Bellows et al, 1981; Ehrlich and Rajaratnam, 1990) in mechanically-loaded environments (Guidry and Grinnell, 1985; Hinz et al, 2001; Mudera, 2000). FPCLs have also been used to investigate the effects of various growth factors on fibroblasts (Grinnell et al, 1999; Schreiber et al, 2001) as well as the behaviour of individual fibroblasts as they contract their environment while moving through the ECM (Roy et al, 1997, 1999). Further details on the range of experiments involving FPCLs, and their clinical utility, are provided in the reviews by Dallon and Ehrlich (2008), Grinnell (2003), and Ehrlich and Moyer (2013).

FPCLs are typically classified according to the mechanical set-up that is used to create them. By this reckoning, there are three main types of FPCLs: free-floating, attached and stress-relaxed.

Free-floating FPCLs were introduced by Bell et al (1979) and are prepared by polymerising a collagen gel with fibroblasts. This could be done either in a bacteriological dish, to which the gel adheres poorly, or in a tissue culture dish, in which case the gel has to be detached after a certain time. The image displayed in Fig. 3, taken from Kelynack (2009), shows an example of a free-floating FPCL. In such lattices, the fibroblasts project a dendritic network of extensions and the tension is distributed isotropically (Grinnell et al, 2003; Rhee and Grinnell, 2007).

Fibroblasts in free-floating lattices can generate significant traction forces (Majno and Joris, 2004; Grinnell, 2000), which reorganise the matrix and lead to contraction. While this reorganisation does not orient collagen fibrils in any particular direction, it can still cause these lattices to contract to as little as a tenth of their initial lateral (or vertical) extent (Bell et al, 1979; Steinberg et al, 1980; Grinnell and Lamke, 1984; Guidry and Grinnell, 1985), even in the absence of protomyofibroblast cells, which can exert greater forces. This contraction gives rise to a mechanically relaxed tissue that resembles dermis, and it has hence been proposed that such lattices can be used to describe the earliest stages of wound healing, before inflammation and tissue stress have activated the differentiation of fibroblasts into myofibroblasts (Grinnell, 1994). Such FPCLs are observed to remain disk-shaped throughout the contraction process, which involves a reduction in thickness as well as diameter, although the edges of the disk are observed to eventually curl up (Bell et al, 1979).

Attached FPCLs are fibroblast-populated lattices that are polymerised in a tissue culture dish, to which the gel attaches firmly. A consequence of this experimental arrangement is that the lattices decrease in thickness but not in lateral area (Grinnell and Lamke, 1984; Guidry and Grinnell, 1985). The tension in such lattices is distributed anisotropically, while fibroblasts develop an elongated bipolar appearance, orienting themselves along the lines of tension (Stopak and Harris, 1982; Bellows et al, 1982; Grinnell, 1994; Tamariz and Grinnell, 2002). This reorganisation causes collagen fibrils to become oriented in the same plane as the substrate, which in turn gives rise to mechanical loading within the matrix. The contraction of such lattices, which involves a reduction in thickness alone, gives rise to a mechanically stressed tissue resembling granulation tissue, and it has therefore been proposed that such lattices can be used to model the early stage of wound healing when the granulation tissue begins to develop and exert stresses on its environment (Grinnell, 1994). The rate and extent of contraction of the lattice is similar to that in experiments where fibroblasts are embedded within such lattices  (Grinnell and Lamke, 1984; Guidry and Grinnell, 1985, 1986). Fibroblasts in such lattices organise a fibronectin matrix and develop prominent actin stress fibres (Farsi and Aubin, 1984; Mochitate et al, 1991; Halliday and Tomasek, 1995), which indicate that some fibroblasts have differentiated into the protomyofibroblast phenotype. It has also been observed that fibroblasts in restrained matrices develop fibronexus junctions (Tomasek et al, 2002), which allow TGF-$\beta$, if present, to further stimulate the differentiation of protomyofibroblasts into the myofibroblast phenotype (Arora et al, 1999; Vaughan et al, 2000).

A Stress-relaxed FPCL is prepared by polymerising collagen lattice, allowing it to attach to a tissue culture dish for a set period of time and then detaching it (Tomasek et al, 1992). In such lattices, tensile stress develops while the matrix is anchored and these stresses are relieved via a sudden smooth muscle-like contraction when the matrix is released, as the cell extensions collapse and the stress fibres disappear (Mochitate et al, 1991; Tomasek et al, 1992; Lin et al, 1997). It has been proposed that stress-relaxed FPCLs can be used to model scar tissue, or the transition from granulation tissue to replacement dermis in wound healing or in tissue repair (Carlson and Longaker, 2004).

As discussed earlier, the contraction of an FPCL results from a fundamental reorganization of the lattice structure, and is hence effectively “permanent”. We now discuss the mechanisms that have been proposed to explain the contraction of FPCLs.

2.2 Mechanical properties

It is known that the degradation and replacement of collagen fibres is not a significant mechanism of contraction in attached lattices. This was demonstrated experimentally using radiolabelled collagen, where it was found that the concentration of proteins in the gel remains constant during the reorganisation process (Guidry and Grinnell, 1985). Indeed, it seems plausible to assume that contraction in all types of FPCLs is largely a result of the rearrangement of pre-existing collagen fibres by fibroblasts. Even still, the precise mechanism of lattice contraction is thought to vary, depending on the density and mechanical state of the lattice.

The prime mechanism of stressed lattice contraction is believed to be cell contraction, which is associated with the protomyofibroblast and myofibroblast phenotypes that are prevalent in such environments (Dallon and Ehrlich, 2008). When these cells contract, they pull on the surrounding lattice, causing it to contract with them. It has been suggested that contraction in free-floating lattices of moderate cell density occurs through a cell traction mechanism (Dallon and Ehrlich, 2008), in which cell locomotion results in the compacting of collagen fibres by bundling thin fibrils (Harris et al, 1980; Ehrlich, 2003). As discussed in Dallon and Ehrlich (2008), the contraction of free-floating lattices of high cell density is believed to occur through a cell elongation and spreading mechanism, in which fibroblasts pull collagen fibrils towards them, thus compacting the gel.

A salient feature of FPCLs is that they provide a testbed for investigating the how the interplay between fibroblasts and the collagen matrix can cause the latter to permanently contract. Under the action of fibroblasts, the density of neighbouring fibrils increases locally and consequently the volume of the collagen lattice decreases (Grinnell, 2003). An example of such behaviour, taken from an experiment by Kelynack (2009), is shown in Fig. 3.

The observation that the lattices only re-expand partially crucially indicates that contraction is not simply the elastic response of the collagen lattice to traction forces applied by cells. Unlike other forms of remodelling, such as the hardening of the human eye lens described by Augusteyn (2010), which could be modelled using a time-dependent stress-strain relationship, FPCL contraction involves some form of cell-induced ‘plastic’ behaviour. That is, the active lattice remodelling by cells causes the unloaded state of the lattice to change over time in a manner analogous to classical plasticity. Hence, a central aim of this paper is to develop a mathematical framework for FPCL contraction that takes into account the evolving unloaded state. Our approach quantitatively captures key features of the contraction process, which would be impossible to achieve with, for example, a Kelvin-Voigt viscoelastic constitutive law. We now briefly outline some of the previous mathematical approaches to modelling the behaviour of FPCLs.

2.3 Modelling FPCL contraction

Most previous models of FPCL contraction have incorporated a viscoelastic framework. The earliest such model was developed by Moon and Tranquillo (1993), who adapted the Tranquillo and Murray (1992) model of dermal wound healing to describe the contraction of a collagen microsphere. This model could not, however, account for the permanence of matrix contraction. An alternative framework, which addresses this issue, but which is only valid for small displacement gradients, was developed by Barocas et al (1995) who replaced the Kelvin-Voigt constitutive law in the Moon-Tranquillo model with a Maxwell constitutive law. In contrast, Ferrenq et al (1997) retained the linear Kelvin-Voigt constitutive law, but restricted their focus to situations in which the displacement gradient is small and linear theory is valid. A different approach, that used the theory of mixtures to describe the interaction between the fibrous lattice and the permeating fluid medium, was the biphasic model of collagen lattice contraction developed by Barocas and coworkers (Barocas and Tranquillo, 1994, 1997). This model takes into account several important effects, including the partial expansion of collagen lattices after cell traction stresses are removed. Subsequently, Barocas and coworkers have used similar approaches to describe several experiments on collagen lattices (see, for example, Chandran and Barocas (2004); Knapp et al (1999); Schreiber et al (2003)), using models that include fibroblast traction as an additive contribution to the total stress. Recent approaches to modelling FPCL contraction include models by Marquez, Zahalak and coworkers (for instance (Marquez et al, 2005; Pryse et al, 2003; Zahalak et al, 2000)), who considered individual cells and used Eshelby’s solution to describe the local strain fields and Green et al (2013), who considered spherical-symmetric collagen lattices and treated the gel as a compressible Stokes fluid.

While these models have significantly advanced the understanding of this process, there have not yet been any models that explicitly take into account the continual mechanical restructuring of such lattices. In the next section we introduce a morphoelastic model, based on the concept of an “effective strain” that can be used to capture the remodelling of FPCLs.

3 A morphoelastic model for the contraction of fibroblast-populated collagen lattices

In this section we develop a mechanical model that can be used to describe the contraction of a collagen lattice by fibroblasts. As described in standard continuum mechanics texts (for example Gonzalez and Stuart (2008)), it is common to work in either a ‘Lagrangian’ (or ‘material’) coordinate system, in which each particle is labelled according to its position in the initial configuration of the body, or in an ‘Eulerian’ (or ‘spatial’) coordinate system, in which each particle is labelled according to its current position. Holding the Lagrangian coordinate constant corresponds to observing a single particle, while holding the Eulerian coordinate constant corresponds to observing a single point in space.

As we discuss in Sec. 3.2, it is often appropriate to treat contracting FPCLs as $1$-D bodies. Indeed, Ferrenq et al (1997) constructed a model of FPCL contraction in $1$-D Cartesian coordinates and various authors have assumed radial symmetry to develop $1$-D descriptions of FPCL contraction (Moon and Tranquillo, 1993; Ramtani et al, 2002; Ramtani, 2004). We hence derive a one-dimensional model, based on the following equation (derived explicitly in Appendix A) that is appropriate for describing the mechanical behaviour of morphoelastic solids with small effective strain

[TABLE]

where $e^{E}$ is the Eulerian strain, $v$ is the velocity of a point in the material and $g(x,\,t)$ describes the rate of growth. As discussed in more detail in Appendix A, Eulerian strain is taken here to be the effective Eulerian strain, a local, dimensionless measure of the difference between the current state and the zero stress state, while the rate of growth is a measure of the rate at which the zero stress state becomes larger over time. In Sec. 3.1, we take advantage of the simplifications that can be made to eq. (2) when the deformation of the lattice might be large, but the difference between the current state and the stress-free state is always small. In particular, we introduce a form for the growth function, $g(x,\,t)$, based on a plausible mechanism for the contraction of the lattice by fibroblasts. This ultimately leads to a full model of FPCL contraction, which we present in Sec. 3.2.

3.1 A cell-based contraction model

In the experiments performed by Guidry and Grinnell (1985), lattices were contracted by fibroblasts and then allowed to undergo a partial re-expansion after reorganisation is inhibited. This re-expansion was faster than the contraction process, but far from instantaneous, indicating that the time scale associated with the viscous relaxation of collagen lattices should not be ignored. It is hence preferable for us to use a Kelvin-Voigt viscoelastic constitutive law to relate stress and effective strain instead of the purely elastic law, namely $\sigma=E\,e^{E}$, described in Appendix A.2.

Additionally, the large deformations associated with lattice contraction may cause significant changes to the elastic properties of the collagen: We expect the collagen to become stiffer as it becomes denser. Following Ramtani et al (2002) and Ramtani (2004), we therefore propose that the elastic modulus of the collagen should be a function of collagen density. Incorporating viscoelasticity and the changing elastic modulus (but ignoring the activity of cells), this means that an appropriate constitutive law for a collagen lattice will take the form

[TABLE]

where $\mathcal{E}(\rho)$ is the elastic modulus, $\rho$ is the collagen density and $\mu$ is the collagen viscosity.

Next, we develop an expression for $g(x,\,t)$ that captures how cells actively rearrange the fibres of the collagen lattice by applying traction – a microscopic mechanism that causes a modification in size of the zero stress state. Our approach is to imagine the collagen lattice as a $1$-D body containing evenly-spaced fibroblasts, each of which effectively acts as a force dipole by pulling on the collagen lying on either side of it and compressing the collagen directly under it. We assume that each fibroblast rearranges the collagen directly under itself, thus evolving the lattice to a permanently compressed state. Thus, $g(x,\,t)$ will be negative and directly proportional to the compressive strain in the lattice under the cells. Specifically, we assume that $g(x,\,t)$ is proportional to the amount that the strain in the region covered by the cells $e_{\text{cells}}$ exceeds (i.e. is more negative than) a critical level of contraction, $-\hat{e}_{\text{crit}}$.

It has been observed that even though fibroblasts lead to substantial deformation of the environment, large strains are confined to regions around each cell (Sander, 2013), which results in heterogeneity in collagen density. Hence, to realistically describe the interplay of cells with their environment, one would need to assume that gel compaction is not spatially homogeneous, along the lines of the models by Evans and Barocas (2009) or Stevenson et al (2010). This approach would introduce a range of complications that do not appear to be physically relevant to FPCLs; as noted in Grinnell and Lamke (1984), apparently homogeneous compaction of a gel is observed even when cells are cultured purely on the gel surface.

For simplicity and brevity, in the analysis that follows we assume that the FPCL consists of a periodic array of identical units, each of which contains a single cell at its centre (see Figure 4), and hence that the fibroblast density is spatially constant. We assume that the traction forces applied by cells create regions of the lattice (under the cells) that are more compressed than other regions of the lattice (outside the cells). However, we also assume that the zero stress state remains homogeneous in space even as it evolves over time. This modelling choice is made in order to account for some effects of cell mobility. We assume that cells make local changes to the zero stress state, and then move nearby and repeat this process. The overall effect of this will be for the changes to the zero stress state to be homogenised through space.

Using our assumption of a periodic array of cells, we construct a local coordinate system in which a given unit extends from $-l_{u}$ to $l_{u}$ and the cell extends from $-l_{c}$ to $l_{c}$. Thus, the cell density, $n$, is proportional to $l_{c}/l_{u}$. We also assume that the zero stress state is uniform across the unit, and that if the unit was stress-free throughout it would extend from $-l_{z}$ to $l_{z}$. Hence, the average effective strain throughout the unit is

[TABLE]

Note that a further advantage of the assumption of uniform cell density across the array is that this expression for strain will be equivalent to that measured in a representative volume element containing many cells. We shall therefore use this measure of strain in our macroscopic model.

As the cell in each unit pulls the lattice in towards itself, the collagen under the cell will be under less tension (or more compression) than that over the rest of the unit. More formally, a cell in a given unit can be treated as a pair of body forces, both of magnitude $\tau_{c}$; at $-l_{c}$ the cell pulls the collagen lattice to the right, while at $l_{c}$ the cell pulls the lattice to the left. Neglecting inertial terms, the momentum balance equation within a single unit is then

[TABLE]

where $\sigma_{u}$ is the local stress and $\delta$ is the Dirac delta function. Integrating this equation with periodic boundary conditions, we find that the stress within the unit is

[TABLE]

where $\sigma_{\text{out}}$ represents the stress in the regions outside the cells, which we assume to be governed by the constitutive law (3), so that

[TABLE]

Since each cell must apply a balanced pair of body forces to the lattice, it follows that either the boundary of the lattice occurs in a region where $\sigma=\sigma_{\text{out}}$, or there is an additional surface traction term that accounts for the effects of the end of the cell. In either case, we find that it is most convenient to express macroscopic boundary conditions in terms of $\sigma_{\text{out}}$: a free-floating lattice will have $\sigma_{\text{out}}=0$ on the boundary, while a spring connected to the boundary (as in, for example, Marenzana et al (2006)) should be expressed as a relationship between $\sigma_{\text{out}}$ and the displacement at the boundary.

From the definition of strain above, we can express $\tau_{c}$ as

[TABLE]

where $\sigma_{\text{cells}}$ is the stress in the region covered by the cell. Assuming that viscous stresses are uniform (or negligible) between $l_{c}$ and $-l_{c}$, and using a constitutive law for $\sigma_{\text{cells}}$, we obtain the expression

[TABLE]

As mentioned earlier, we assume that the growth rate $g(x,\,t)$ is dependent on the physical contraction experienced by the lattice directly under the cells and is proportional to $-(e_{\text{cells}}-\hat{e}_{\text{crit}})$. Moreover, the rate of contraction is proportional to the cell density, $n$, as this is representative of the proportion of the lattice that is accessible to the fibroblasts. Thus, we obtain the following constitutive law for $g(x,\,t)$:

[TABLE]

where $\theta$ is a constant of proportionality with dimensions of cell density*-1* time*-1* and the positive part operator $(X)^{+}=X$ for any positive $X$ and is zero otherwise. As the condition

[TABLE]

is always satisfied in most practical situations, we neglect the positive part operator in (7) for the rest of our analysis.

Next, we note that since the zero stress state is uniform throughout the unit,

[TABLE]

Combining this with (4) yields the expression

[TABLE]

Using (6), and defining $\sigma_{c}$ as rescaling of $\tau_{c}$ so that $n\,\sigma_{c}=l_{c}/l_{u}\,\tau_{c}$, we obtain

[TABLE]

Substituting into (5), gives us our constitutive relationship between stress and strain:

[TABLE]

Note that the $n\,\sigma_{c}$ term is analogous to the cell traction stress term found in other models of lattice contraction and dermal wound healing (Moon and Tranquillo, 1993; Tranquillo and Murray, 1992; Ferrenq et al, 1997; Tracqui et al, 1995). Indeed, by making $\sigma_{c}$ a function of the lattice density and/or cell density (to incorporate the effects of crowding on traction stress, for example) we can recover identical expressions to those used in these earlier papers. However, our approach differs in that traction stress emerges naturally from the assumption that cells apply body forces to the collagen lattice, rather than by requiring that this stress term be incorporated into the constitutive law. Nevertheless, the consistency between these two treatments of cell traction gives us further confidence that both approaches are valid and useful.

Dropping the subscripts on $e_{\text{avg}}$ and $\sigma_{\text{out}}$, we obtain the following set of equations to describe the mechanical behaviour of a contracting viscoelastic lattice:

[TABLE]

where $n_{l}=n\,l_{u}/l_{c}$ is the constant of proportionality that relates the cell density, $n$, to the cell length per unit length, $l_{c}/l_{u}$.

A full description of the mechanics of a contracting lattice will also require an initial condition on $e$ and two boundary conditions, either on $\sigma$ or $v$, that specify whether the ends of the lattice are tethered or stress-free. We now combine (9) with some assumptions about fibroblast motion and interactions in FPCLs to construct a full morphoelastic model of FPCL contraction.

3.2 The morphoelastic model of FPCL contraction

In the following we assume that the contraction of the lattice is due to the rearrangement of collagen fibrils by fibroblasts alone, with the expectation that the results obtained using this assumption will be quantitatively similar to those obtained using a model that takes into account the rôle of protomyofibroblasts. Additionally, we assume that TGF-$\beta$ is not present in the collagen lattice. As protomyofibroblasts require the presence of TGF-$\beta$ to differentiate into myofibroblasts (Tomasek et al, 2002), we hence do not consider the activity of the latter. Furthermore, while it has been observed that fibroblasts in free-floating lattices become quiescent (Rosenfeldt and Grinnell, 2000) and a small fraction of these cells undergo apoptosis (cell death) (Fluck et al, 1998; Grinnell et al, 1999), we shall neglect these processes in our model.

Experimental data for the contraction of FPCLs are typically obtained by measuring the changing diameter (in the case of free-floating FPCLs) or thickness (in the case of attached FPCLs) of the lattice. The behaviour of such lattices can hence be modelled by considering one spatial dimension, that represents either the radius or thickness of the lattice depending on the type of FPCL being considered. It should be noted that the mechanical model presented in Appendix A assumes a single Cartesian spatial dimension and cannot be easily modified to obtain equations for the mechanics of radially symmetric lattice contraction. As a result, our model is best suited to attached, rather than free-floating, FPCLs. Nevertheless, as seen later in this section, our model is equally successful at fitting data obtained from both these lattice types.

Following the conventions used in Appendix A.1, we use $x(X,\,t)$ to represent the position at time $t$ of a particle initially located at $X$ and $X(x,\,t)$ to represent the initial position of a particle at position $x$ at time $t$. Since the lattice changes in size over time, the domain of interest can be expressed as $0\leq x\leq l(t)$ or, equivalently, $0\leq X\leq l_{0}$, where $l_{0}=l(0)$. As illustrated in Fig. 5, the point $x=0$ is taken to represent either the centre of a floating lattice or the fixed tethering point of an attached lattice. In either case, this point is fixed and we hence find that $x(0,\,t)=X(0,\,t)=0$.

We next define the Eulerian displacement, $u(x,\,t)=x-X(x,\,t)$, and the Eulerian displacement gradient, $w(x,\,t)={\partial u}/{\partial x}$. Following the definitions of $x$ and $X$, it follows that $u(x,\,0)=w(x,\,0)=0$. Moreover, the fact that the centre or tethering point is fixed implies that $u(0,\,t)=0$. Another important kinematic variable is the velocity, which is defined as

[TABLE]

with $v(0,\,t)=0$ at the fixed centre/tethering point.

Next, we introduce expressions for the densities of the ECM, $\rho(x,\,t)$, and the fibroblast cells, $n(x,\,t)$. Since Guidry and Grinnell (1985) observed that collagen synthesis and degradation are insignificant in FPCLs, we propose that contraction proceeds only by rearrangement of the collagen lattice. In other words, we assume that the net creation of each species is zero during the timescale of the experiment. Furthermore, for this study, we assume that there is no migration of fibroblasts. As both species are subjected to passive advection, this assumption allows us to use simple continuity equations for each of the species, which can be solved explicitly (Clement, 1978) to give

[TABLE]

where $n_{0}$ and $\rho_{0}$ are the (spatially uniform) initial densities of fibroblast cells and the ECM, respectively.

We assume that the lattice is initially relaxed ($e(x,\,0)=0$), and that lattice stiffness increases with density, according to the simple linear relationship

[TABLE]

where $\mathcal{E}_{0}$ is the elastic modulus when $\rho=\rho_{0}$ and $k$ is a positive constant. It should be noted that other forms have been proposed to describe the dependence of the elastic modulus on density. For example, Ramtani and coworkers describe this relationship as a power law (Ramtani et al, 2002; Ramtani, 2004). In the absence of further data, however, we use a simple, linear law.

The case of a free-floating FPCL (Fig. 5a) or an attached FPCL (Fig. 5b) can be described by using the stress-free boundary condition

[TABLE]

and the displacement at $x=l(t)$ can be obtained by noting that

[TABLE]

As we shall see, $w(x,\,t)$ is spatially homogeneous in the problem that we analyse. In this situation, it follows that

[TABLE]

Taking the spatial derivative of (10) and rearranging, we obtain the following evolution equation for $w(x,\,t)$:

[TABLE]

In conjunction with the equations for the evolution of the effective strain (9a), the constitutive law for the total stress (9b) and the force balance equation (9c), equations (11)-(15) comprise an Eulerian model for the contraction of collagen by fibroblasts. This system can be non-dimensionalised by introducing the dimensionless variables

[TABLE]

and the dimensionless constants

[TABLE]

where, for ease of comparison with experiment (see later), we choose $h=1$ hr.

3.3 Lagrangian description of the model

As this Eulerian model has a moving boundary at $x=l(t)$, we now adopt a Lagrangian coordinate system to obtain a model for a system in which the right hand boundary is fixed. Conveniently, we find that transforming our equations to Lagrangian variables is equivalent to converting to characteristic variables: we ultimately obtain a system where the only differentiation is with respect to time. In the derivation that follows, it is important to distinguish between partial time derivatives with the Eulerian spatial coordinate held fixed and partial time derivatives with the Lagrangian spatial coordinate $X$ held fixed. Hence, we introduce the ‘Lagrangian time variable’, $T$, and we transform the entire problem from $(x,\,t)$ coordinates to $(X,\,T)$ coordinates.

Firstly, we introduce Lagrangian displacement $U$ and Lagrangian velocity $V$, which are defined as

[TABLE]

These are both identical to $u(x,\,t)$ and $v(x,\,t)$ with just a transformation of the independent variables (see Appendix A.2 for details). Now, we define the Lagrangian displacement gradient, $W(X,\,T)={\partial U}/{\partial X}$, which is different from $w$ since differentiation with respect to $x$ is different from that with respect to $X$. As described in Appendix B, however, the definitions of $W$ and $V$ lead to the identities

[TABLE]

and hence $w=W/(1+W)$. Based on these identities, it follows that

[TABLE]

where $\Phi=\phi\,(1+W)$ and where $\phi$ is a general scalar quantity. This motivates the introduction of Lagrangian variables $N=n\,(1+W)$, $R=\rho\,(1+W)$, $S=\sigma\,(1+W)$ and $E=e\,(1+W)$, so that all of the advective time derivatives simply become partial time derivatives. From (11) and (12), we obtain the following exact expressions for the scaled densities $N$ and $R$ in the ‘normal’ contracting situation

[TABLE]

Note that the situation in which lattice reorganisation is inhibited by the addition of cytochalasin D at some time, $T_{\rm inh}$, can be described by taking

[TABLE]

From (15), (9a), (9b) and (9c) we obtain the equations

[TABLE]

where $\mathscr{E}(W)$ is

[TABLE]

Using the stress-free boundary condition (13), it follows that

[TABLE]

As $W$ is independent of $X$, it follows that $w$ is independent of $x$ and hence (14) is valid. Moreover, the fact that $W$ and $E$ are both independent of $X$ means that we can rewrite our system of PDEs as a pair of coupled ODEs for the displacement gradient and strain:

[TABLE]

Although this model has five dimensionless free parameters: $\bar{\mu}$, $\bar{k}$, $\bar{\tau}$, $\bar{e}_{\rm{crit}}$ and $\bar{\theta}$, as we show in the next section, we can fix three of these parameters using heuristic arguments, leaving us with a two-component model with two free parameters. Finally, on non-dimensionalising (14) and converting to Lagrangian variables, we find

[TABLE]

Thus, our transformed ODE model gives easy access to the evolving length of the lattice, the physical variable variable most easily observed in experiments.

3.4 Comparison with experimental data

In the following, we numerically integrate the system (22)-(23) and compare the results with previously obtained data. In particular, we consider the experiments of Bell et al (1979), Talas et al (1997) and Feng et al (2003), which used free-floating FPCLs, as well of Guidry and Grinnell (1985), which was performed using an attached FPCL and in which lattice reorganisation was halted at various times. We first estimate values for the parameters $\mathcal{E}_{0}$, $F_{\text{cell}}$ and $\sigma_{0}$, which we then use to fix the parameters $\bar{\mu}$, $\bar{\tau}$ and $\bar{e}_{\rm crit}$ for each case (values listed in Table 1).

The initial stiffness of a collagen gel $\mathcal{E}_{0}$ was measured by Knapp et al (1997, 1999) to be $1.185$ kPa. However the precise value of this quantity has been found to vary widely, depending on the extent of cross-linking (Discher et al, 2005). As the initial stiffness of the gel was not measured in any of the experiments under consideration, we assume for the purposes of our simulations that the gel used in each case was a lightly cross-linked collagen lattice, and use the approximate value of the stiffness of such lattices (Discher et al, 2005), namely $\mathcal{E}_{0}=1$ kPa.

Next we note that the total cell traction force will be a product of the applied contraction force per cell $F_{\text{cell}}$ and the number of cells in the lattice, $C$. On dividing this quantity by the cross-sectional area of the lattice, $A$, we obtain an expression for initial cell traction stress:

[TABLE]

The values of $C$ and the initial lattice diameter $d$ (from which we can calculate $A$) for each of the four experiments are listed in Table 1. Although the value of $F_{\text{cell}}$ was not measured in any of the four experiments, there have been several other experiments in which the contraction forces exerted by the cells have been determined using cell-populated lattices attached to force monitors, and a wide range of values of $F_{\text{cell}}$ has been reported, ranging from from $0.1$nN/cell (Eastwood et al, 1996) to $1000$nN/cell (Wakatsuki et al, 2000). Here, we will use the result from Kolodney and Wysolmerski (1992), $F_{\rm cell}=500n$N/cell, which was obtained using fibroblast-populated matrices similar to those used in the experiments by Bell et al (1979), Talas et al (1997), Feng et al (2003) and Guidry and Grinnell (1985).

Now, we note that equation (9b) implies that the initial cell traction stress can also be expressed as

[TABLE]

Rearranging the equation for $\bar{\tau}$ in (16) and using (25) with this alternative expression for $\sigma_{0}$, we find that

[TABLE]

Thus, the value of $\bar{\tau}$ can be estimated for each experiment, given the values of $C$ and $A$.

We next note from (16) that the viscosity of the lattice is given by

[TABLE]

For the purpose of our simulations, we assume that $\mu=1.0\times 10^{7}$Pa sec, a value in agreement with Knapp et al (1997) who measured the shear viscosity of a collagen gel to be $1.24\times 10^{7}$Pa sec. As we assume that the initial gel stiffness is $\mathcal{E}_{0}=1$kPa, we hence set $\bar{\mu}=2.78$ for all experiments.

Finally, we note that although $\bar{e}_{\rm crit}$ is difficult to determine experimentally, it can be observed from (16) that it is related to the cell density. Hence, we make the heuristic assumption $\bar{e}_{\rm crit}=C/(C_{f}\,A)$, where $C_{f}$ represents the characteristic number of fibroblasts in an FPCL. In the following, we choose $C_{f}=2\times 10^{6}$, which is a reasonable estimate that lies within the range of values of $C$ for the experiments considered here. It is interesting to note that the above assumptions yield a simple relationship $\bar{e}_{\rm crit}=10^{-3}\,\bar{\tau}$.

The remaining parameters $\bar{\theta}$ and $\bar{k}$ are difficult to estimate from experiments. In particular, $\bar{\theta}$ arises from the constitutive law for $g(x,t)$ and there are no constraints on its value. Hence, in the following we obtain and use “best-fit” values of $\bar{\theta}$. Moreover, we note the power law exponent for the relationship between elastic modulus and density used Ramtani and coworkers used (Ramtani et al, 2002; Ramtani, 2004) is, to a first order approximation, identical to $\bar{k}$. As this exponent was taken to be a free parameter in these previous approaches, we obtain and use best-fit values of $\bar{k}$.

We now qualitatively describe the behaviour of the FPCL in each of the four experiments under consideration by using the values displayed in Table 1 and varying the free parameters $\bar{\theta}$ and $\bar{k}$ using the MATLAB subroutine fminsearch, such that we minimise the difference between the experimental values for the fraction of the original length (diameter or thickness) and the values of $l(t)$, given by (24) at the corresponding measurement times. Specifically, we find the values of $\bar{\theta}$ and $\bar{k}$ that yield the lowest values of the quantity $\chi^{2}$, defined as:

[TABLE]

where $i$ correspond to all the experimental data points being considered, $\mathscr{E}_{i}$ are the values of $l(t)$ for these data points and $\mathscr{O}_{i}$ are the corresponding values of $1+W(T)$ obtained from the simulations. As shown in Table 2, we find that the $\chi^{2}$ values, obtained by using the optimal values of $\bar{\theta}$ and $\bar{k}$ for each case, are in the range $0.004-0.059$, indicating a good fit between theory and experiments.

To test the robustness of these fits, we determine the values of $\chi^{2}$ obtained by simulating our model for a wide range of choices of $\bar{\theta}$ and $\bar{k}$. In Fig. 6, we display contour plots that indicate how $\chi^{2}$ varies with $\bar{\theta}$ and $\bar{k}$ for each of the four experiments mentioned earlier. In addition, for each case we display the location of the “optimal” choice of $\bar{\theta}$ and $\bar{k}$, obtained through the minimization procedure described above. It can be observed from this figure that small perturbations (around $\pm 10\%$) in the value of either of the two parameters from the global minima will not substantially change the resulting values of $\chi^{2}$. Indeed, as seen in Table 2, the corresponding $\chi^{2}$ values in these cases are of a similar order of magnitude to those obtained when using the best-fit values of $\bar{\theta}$ and $\bar{k}$. This implies that the global minima displayed in Fig. 6 can be used to obtain robust fits of the simulations to each of the four datasets.

It is important to note that our choice of the global minima of $\chi^{2}$ is in order to provide a unifying criterion for the selection of $\bar{\theta}$ and $\bar{k}$. Although the best-fit values listed in Table 2 vary over a wide range, there is a large region of the ($\bar{\theta}$,$\bar{k}$)-space over which we could obtain fits that are essentially as good (see Fig. 6). Thus, one could have, in principle, chosen very similar values of $\bar{\theta}$ and $\bar{k}$ for each of the three free-floating FPCL experiments and have still obtained good fits in each case. While such a choice might potentially be more realistic, as one would not expect much variance of these parameters across different free-floating FPCLs, it is difficult to justify the use of any pair of ($\bar{\theta}$,$\bar{k}$) over another. To this end, we restrict our attention to the best-fit values in Table 2.

We next examine in detail the best-fits to the individual datasets. The data for the change in the fraction of the original lattice extent for four different experiments are shown in Fig. 7. It is important to note from Table 1 that the lattices for the four chosen experiments have very different spatial extents, which explains the differing time scales for contraction. In fact, it appears that the rate at which the fraction of the original extent of a lattice shrinks is inversely proportional to its original extent. This suggests that the contraction rate of the lattices are related to their initial sizes.

We begin by considering the experiments performed by Bell et al (1979) that yielded the first reported observation of FPCL contraction. In these experiments, different numbers of fibroblast cells were embedded in free-floating collagen lattices of varying ECM density, which contained fetal bovine serum. The fraction of the original length was then measured at various times over the course of several days. In the following, we consider their set of results that were obtained using $7.5\times 10^{6}$ fibroblast cells embedded in a lattice of initial diameter $53$mm. It was observed that this lattice subsequently contracted to around $10\%$ of its original diameter in around $9$ days. We find that the best fit of our numerical results to the experimental data can be obtained by using $\bar{\theta}=0.89$, $\bar{k}=106.62$ (see Fig. 7a).

We next consider the experiments of Talas et al (1997), in which the difference between the effects of normal and recessive dystrophic epidermolysis bullosa fibroblasts were investigated. For the case of normal fibroblasts, $2.5\times 10^{5}$ cells were embedded in a lattice with a diameter of approximately $2.18$ cm. It was observed that the lattice contracted to around $37\%$ of its original area in around $7$ days. We find that the best fit of our numerical results to the experimental data can be obtained by using $\bar{\theta}=5.10$, $\bar{k}=1294.92$ (see Fig. 7b).

Next, we consider the experiments of Feng et al (2003), in which the mechanical properties of contracted collagen gels were investigated. Here, $2.5\times 10^{6}$ fibroblast cells were embedded within a lattice of initial diameter $100$mm. It was observed that the lattice contracted to around $13\%$ of its original diameter by the end of the experiment, with most of the contraction occurring within the first few days. In this case, we find that the best fit of our numerical results to the experimental data can be obtained by using $\bar{\theta}=16.14$, $\bar{k}=183.85$ (see Fig. 7c).

Finally, we consider the results of Guidry and Grinnell (1985), obtained when $10^{5}$ fibroblasts were placed on an attached lattice of diameter $12$mm. It was observed that this lattice contracted to about $24\%$ of its original thickness in around $1$ day. As the presence of protomyofibroblasts was not reported in these experiments, it is likely that the observed contraction is primarily due to the activity of fibroblasts and hence our model can be used to approximate this behaviour. We simulate the effect of adding cytochalasin D by allowing $N(X,\,T)$ to take the form (20). Results obtained using $\bar{\theta}=0.33$, $\bar{k}=2.97$, which gave the best fit to the experimental data, are shown in Fig. 7d. We find that in addition to the contraction of this gel, our model can capture the observed partial re-expansion.

4 Discussion

In this work, we develop a $1$-D morphoelastic model to describe the evolution of the diameter (width) of a free-floating (attached) collagen lattice that is contracted by fibroblasts. We fit numerical solutions of the model to previously obtained experimental data by optimizing two free parameters. In addition to being able to capture the contraction of the lattice, our model closely describes the partial re-expansion of the lattice observed in the experiment by Guidry and Grinnell (1985). Hence, this is to our knowledge the first mathematical model that explains how the permanent contraction of such lattices partially persists when reorganisation is inhibited by killing the fibroblasts or by otherwise preventing them from altering the lattice. We achieve this by explicitly taking into account the continual changes to the zero stress state of a material in response to a prescribed rate of growth and explicitly describes the evolution of this state. The partial re-expansion of the lattice is an example of a class of phenomena related to changes in the underlying tissue structure. As the theoretical framework developed here is flexible and versatile, it could potentially be used to model a range of other biological processes that involve internal remodelling of a tissue’s mechanical structure, by making certain assumptions about the dependence of growth on other physical parameters.

Note that our framework for stress and strain is developed in Eulerian coordinates (in Appendix A). Although a Lagrangian framework can lead to equations with fixed boundaries and a useful variational structure (see, for example Roberts (1994) and Gonzalez and Stuart (2008)) it is typically only appropriate in situations where the zero stress state does not evolve. However, as described in Yavari (2010), the introduction of a changing zero stress state leads to new terms that need to be carefully incorporated into the energy balance equation. Moreover, in cases where the zero stress state is close to the current state, but quite different from the initial state, Eulerian coordinates enable aspects of small deformation theory, such as the linear stress-strain relationship, to be applied.

On deriving the relations for stress and strain in Eulerian coordinates, we then revert to Lagrangian coordinates (in 3.2) in order to avoid the problems inherent in a system with a moving boundary. To do this, we exploit the fact that a conversion to Lagrangian coordinates corresponds to a conversion to characteristic variables (as diffusion is absent in our system), and we are thus able to simplify our model to a system of ordinary differential equations.

It should be noted that we have assumed in our model that the zero stress state is always uniform within each unit and that the expression for $g(x,\,t)$ also gives a uniform rate of change to the zero stress state. Against this, it might be argued that the only region where the zero stress state is changing (due to the direct effect of fibroblast action) is the region where $-l_{c}<x<l_{c}$. However, experiments have been performed where cells were cultured on top of lattices instead of throughout lattices, and these still showed relatively uniform contraction (Grinnell and Lamke, 1984; Guidry and Grinnell, 1985, 1986). This indicates that while cell migration may be slow (Grinnell and Lamke, 1984; Ascione et al, 2016), there is sufficient cell movement (or perhaps cells rearrange collagen over sufficiently large volumes) for it to be reasonable to assume that the change in the zero stress state is uniform in space.

Despite the fact that cell movement is quite limited over the time scale of FPCL contraction, a natural extension to our model would be to consider the case where cell density is nonuniform, and where there is active cell movement over the time-scale of contraction. Extensions like this would be necessary for developing models of wound healing based on the principles presented in this model. When cell motility is taken into consideration, the transformation to Lagrangian coordinates would no longer convert the partial differential equations of the original model to ordinary differential equations, and it is likely that it would be more straight forward to solve the original Eulerian equations.

The morphoelastic approach detailed in this paper, where an evolution law is developed for the effective strain but a conventional stress-strain relationship is used, could potentially be applied to several other biological processes including soft tissue growth, arterial remodelling, aneurysm development, morphogenesis, initiation of stretch marks and solid tumour growth. Furthermore, the permanent contraction observed in some pathological scars instead suggests that morphoelastic changes to the underlying extracellular matrix are particularly significant during the wound healing process. Hence, an important potential application of this theory is the development of a mechanochemical model of dermal wound healing. Indeed, it would be fruitful to build on recent theoretical models of the wound healing process that have considered aspects of the interplay between fibroblasts, keratinocytes and growth factors (Menon et al, 2012), the complementary rôles of TGF-$\beta$ and tissue tension (Murphy et al, 2011a, 2012), the interplay between these components and ECM deformation (Valero et al, 2014), and their rôle in hypertrophic scar formation (Koppenol et al, 2017b) and wound retraction (Koppenol et al, 2017a).

In addition, there are several possible ways in which this model could be extended to describe other observed aspects of the behaviour of FPCLs. For instance, there have been experiments in which force measurement devices were attached to the collagen lattice (Kolodney and Wysolmerski, 1992; Brown et al, 1996; Marenzana et al, 2006) in order to measure the total contractile force exerted by the cells. In each case, the force measurement device provides a finite spring-like resistance to the contraction of the lattice, and so the contraction process could be described using our model by replacing the stress-free boundary condition in (13) with one that relates the stress at $x=l(t)$ to the displacement at the contracting edge. However, as the increased elastic tension in such a lattice could lead fibroblasts to modulate into protomyofibroblasts (Tomasek et al, 2002) that increase the stress, it may also be necessary to make the cell-associated stress, $\sigma_{c}$, a function of elastic stress. Indeed, this modulation can also occur in stress-relaxed FPCLs (Tomasek et al, 1992), which causes the lattice to contract rapidly when released. In the additional presence of TGF-$\beta$, protomyofibroblast cells in such lattices will further modulate into myofibroblasts (Desmoulière et al, 1993; Tomasek et al, 2002; Gabbiani, 2003; Desmoulière et al, 2005). In order to capture the effect of these highly contractile cells, our model could be modified via the addition of a new species for protomyofibroblasts (or myofibroblasts), a new chemical species, TGF-$\beta$, and including a conversion term in the equation for the fibroblasts. Furthermore, as mentioned earlier, the effect of heterogeneity in gel compaction could also be investigated using the theoretical framework that we describe.

However, it is important to note that although our morphoelastic model closely captures the important features of the permanent contraction of FPCLs, the approach used here can only truly be justified in $1$-D Cartesian co-ordinates (see Appendix A.1 for more details). Consequently, although our model can provide a good description of attached lattices, a complete description of cylindrical, free-floating lattices would require significant modifications to the constitutive laws used for the mechanics of the lattice (albeit with only minor adjustments to the reaction-diffusion equations). This could still be achieved using the same principles described in Sec. 3.1 if appropriate modifications are made; for example, we could consider a three-dimensional array of units, each containing a single cell, and each cell could be treated as a sphere of body stresses pulling inward. By homogenising such a system, we could develop a three-dimensional constitutive law for growth analogous to (7); however, any such rule would be significantly more complicated than the one-dimensional rule developed in this manuscript.

An additional complication in this regard is that there remains considerable ambiguity surrounding the precise specification of $3$-D laws for the evolution of the zero stress state (for an overview on recent advances in this area, see Kuhl (2014)). As noted in Ambrosi et al (2011), there has been little success in using thermodynamic arguments to develop general frameworks for morphoelasticity. Furthermore, there are uniqueness issues surrounding the multiplicative decomposition of the deformation gradient and it is difficult to ensure that any phenomenological evolution law is appropriately observer-independent. Although some effort has gone into resolving these problems (especially in the engineering literature – see, for example, Lubarda (2001) and Xiao et al (2006)), the resulting models are often densely expressed and difficult to apply to biological morphoelasticity. As noted in chapters $4$ and $5$ of Hall (2008), some progress can be made by considering possible three-dimensional generalisations of (31), but there is a need for further work in this area. We also note that there is still interest in developing constitutive models that could provide better mathematical descriptions of remodelling in soft biological tissues (see, for example, Comellas et al (2016)).

Despite these challenges, the morphoelastic framework presented in this manuscript has significant advantages over other approaches to describing biological remodelling. In particular, phenomena like the permanent contraction of a collagen lattice, which our model can describe, are inaccessible to classical Kelvin-Voigt models and are very different from the stress-induced contraction observed in Maxwell models.

The central achievement of our work lies in developing a theoretical framework for the mechanics of tissue remodelling that explicitly accounts for the action of cells to change the fundamental structure of a collagen lattice. By defining strain in terms of the difference between the current state and the zero stress state and considering how this strain changes according to both physical deformation and the action of cells, we find that morphoelasticity can be easily incorporated into a model of FPCL contraction. Indeed, this approach allows us to use conventional constitutive laws to relate the stress to the strain. The $1$-D morphoelastic framework described in this paper provides us with a simple and meaningful technique to the describe some of the complexities of biological remodelling, and it has the potential to be extended and used in a wide range of other areas.

Acknowledgements.

This research was primarily supported by the Australian Research Council’s Discovery Projects funding scheme (project number DP0878011). SNM is supported by the IMSc Complex Systems Project ($12^{\rm th}$ Plan). CLH acknowledges the support of the Mathematics Applications Consortium for Science and Industry, funded by the Science Foundation Ireland grant investigator award 12/IA/1683.

Appendix A Effective strain and contraction in a $1$-D morphoelastic body

A.1 The multiplicative decomposition of the deformation gradient

The idea that an observed deformation gradient tensor $\boldsymbol{F}$ could be expressed in terms of elastic and plastic tensor fields through a multiplicative decomposition was first introduced by Bilby et al (1957) and Kröner (1958, 1959). This was further developed by Stojanović et al (1964, 1970) in the context of thermoelasticity and Lee (1969) for the description of metal plasticity at large deformations, while Rodriguez et al (1994) and Cook (1995) were the first to utilize this in the context of biomechanics. The pioneering work of Rodriguez et al (1994) and Cook (1995) was later extended and expanded by Hoger and coworkers (see especially Chen and Hoger (2000)), Goriely and coworkers (Ben Amar and Goriely, 2005; Goriely and Ben Amar, 2007; Goriely et al, 2008), Ambrosi and coworkers (Ambrosi and Mollica, 2004; Ambrosi and Guana, 2007; Ambrosi and Guillou, 2007) and by Vandiver (2009). This biological work has developed alongside applications to thermoelasticity and plasticity, and achievements in these areas have informed each other. This is exemplified by the cross-disciplinary work of Lubarda (2001, 2004) and Rajagopal and coworkers (see, for example, Rajagopal and Srinivasa (2004)). For a comprehensive review of the history and applications of the multiplicative decomposition of deformation gradient, see Lubarda (2004).

In a recent comprehensive review of current work on modelling growth and remodelling, Ambrosi et al (2011) describe a number of applications of the multiplicative decomposition of the deformation gradient, ranging from the remodelling of heart muscle to morphogenesis. The multiplicative decomposition of the deformation gradient is now being used in models of biomechanical phenomena ranging from tissue growth (Ben Amar and Goriely, 2005) to the operation of the heart (Göktepe and Kuhl, 2010; Rausch et al, 2011), although it is important to note that this approach will only be valid when the tissue behaves elastically on the timescale of remodelling (Wyn Jones and Chapman, 2012), and that some authors have general reservations about the use of the multiplicative decomposition on theoretical grounds (Xiao et al, 2006). Indeed, Ambrosi et al (2011) note that there are problems and ambiguities to be resolved when developing appropriate laws to describe the evolution of the growth part of the deformation gradient in response to remodelling. In the context of the present work, it is important to note that many of these difficulties are avoided as we restrict our analysis to the one-dimensional ($1$-D) Cartesian case, although it is still necessary to ensure that the constitutive relation is appropriate for the type of remodelling under consideration.

An accessible introduction to the use of the multiplicative decomposition in biological applications can be found in Goriely and Moulton (2011), which begins with an analysis of a $1$-D growing material that is relevant to the research presented here. It is important to note that a $1$-D body can never be residually stressed: it is impossible to encounter the situation in which the zero stress state cannot be achieved without introducing cuts. Moreover, ensuring observer-independence of time derivatives is much simpler along a single dimension, as it excludes the possibility of a rotating observer. We now develop a simple mathematical framework for modelling the growth or contraction of a $1$-D Cartesian morphoelastic body.

A.2 Strain evolution

The deformation gradient is given by the scalar function

[TABLE]

Following the $1$-D version of equation (1) in Goriely and Moulton (2011), we express this as the product

[TABLE]

where the elastic stretch $\alpha$ is the local size ratio between the current state and the zero stress state, and the growth stretch $\gamma$ is the local size ratio between the zero stress state and the initial state (the growth stretch).

In Goriely and Moulton (2011), the constitutive relation for $1$-D growth is assumed to depend on the rate of growth $g(x,\,t)$:

[TABLE]

While (30) is useful at small deformations (i.e. when $F\approx 1$), it leads to inconsistencies if the current state is significantly different from the initial state, as in the case of FPCL contraction. In order to obtain an equivalent of (30) for large deformations, we first note that $g(x,\,t)$ should satisfy

[TABLE]

That is, $g(x,\,t)$ should be defined with reference to the current configuration, but it should measure the rate of change of the zero stress state of any collection of material particles. It follows from (31) that $g(x,\,t)$ is related to the material derivative of $\gamma(X,\,t)$:

[TABLE]

Note that this reduces to (30) when $F\equiv 1$.

Now, we expect that the stress at any point in the body will be related to the difference between the zero stress state and the current state. A plausible constitutive law that relates the stress, $\sigma$, to the elastic stretch, $\alpha$ is

[TABLE]

where $E$ is the Young’s modulus. This is analogous to Hooke’s law for a linear elastic material, but uses an Eulerian rather than a pseudo-Lagrangian measure of strain, since

[TABLE]

where $\Delta x$ and $\Delta z$ relate to the changes in the current and zero stress states, respectively. In cases where the current state is close to the zero stress state, and hence $\alpha\approx 1$, we see that $e^{E}\approx e^{L}$ and (33) is equivalent to other plausible constitutive laws, such as those in Goriely and Moulton (2011):

[TABLE]

Experimental observations indicate that most of the change in size of a contracting FPCL is due to the permanent rearrangement of fibres by fibroblasts (Guidry and Grinnell, 1985). Hence, it is appropriate to use (33) and assume a linear relationship between stress and strain, rather than the nonlinear model (35). Moreover, (33) has an interesting advantage over (34) and (35), namely that the evolution of Eulerian strain in response to growth can neatly be expressed as an advection equation with a source term that is independent of $e^{E}$.

In order to see this, we substitute (29) into (32) to obtain

[TABLE]

Since $F^{-1}\,DF/Dt=\partial v/\partial x$ where $v$ is the velocity and where $\partial v/\partial x$ is the velocity gradient, it follows that $\alpha^{-1}$ satisfies the equation

[TABLE]

and using $e^{E}\equiv 1-\alpha^{-1}$ we thus obtain the mechanical model for a morphoelastic solid with small effective strain

[TABLE]

Note that thus far we have made the assumption that the relation between stress and strain is purely elastic

[TABLE]

However, as discussed in Sec. 3.1, this can easily be extended to viscoelastic bodies through the use of a Kelvin-Voigt viscoelastic constitutive law. In Sec. 3 of the main text, we use this formulation together with (37) to derive a set of governing equations for the contraction of a $1$-D morphoelastic body.

Appendix B Derivation of the spatial and temporal transformations between coordinate systems

In one spatial dimension, we use $X$ to represent the Lagrangian coordinate and $x$ to represent the Eulerian coordinate. At any given time $t$, there will be a one-to-one mapping from the initial configuration to the current configuration. Thus, we can always write $X=X(x,\,t)$ and $x=x(X,\,t)$; moreover, the fact that particles are not permitted to move through each other implies that ${\partial x}/{\partial X}>0$.

Now, the Eulerian displacement gradient is the spatial derivative of $u(x,\,t)=x-X(x,\,t)$, i.e.

[TABLE]

Similarly, the Lagrangian displacement gradient is

[TABLE]

We thus have the relation

[TABLE]

Using the chain rule, the Eulerian spatial derivative is

[TABLE]

The derivative $\partial T/\partial x$ is equal to zero, and using (38) we obtain the expression

[TABLE]

Furthermore, using (38) and (39), we obtain the following transformation for the spatial derivative

[TABLE]

We similarly use the chain rule to obtain the following expression for the Eulerian temporal derivative

[TABLE]

The derivative $\partial T/\partial t$ is equal to one, and so using (40) we have

[TABLE]

Using (15), we have

[TABLE]

which, on using (39) and the fact that $v(0,\,t)=0$, yields

[TABLE]

Now, from (10), we have $v\,(1-w)=\partial{u}/\partial{t}$. Since $u\equiv U$, using the above expression in conjunction with (39), this yields

[TABLE]

where we have used the definitions of $V$ and $W$. This implies that $v\equiv V$, and so we obtain the following transformation for the temporal derivative

[TABLE]

Now, from (15) we have $\partial{w}/\partial{t}+\partial{(w\,v)}/\partial{x}=\partial{v}/\partial{x}$. Using (39), (41) and (42), this reduces to the following relation between the Lagrangian displacement gradient and velocity:

[TABLE]

We use this expression in Sec. 3.2 of the main text to derive our morphoelastic model of FPCL contraction.