A mathematical model for cell-induced gel contraction incorporating osmotic effects

Biological tissues are composed of cells surrounded by the extracellular matrix (ECM). The ECM can be thought of as a fibrous polymer network, acting as a natural scaffolding to provide mechanical support to the cells. Reciprocal mechanical and chemical interactions between the cells and the ECM are crucial in regulating the development of tissues and maintaining their functionality. Hence, to maintain in vivo-like behaviour when cells are cultured in vitro, they are often seeded in a gel, which aims to mimic the ECM. In this paper, we present a mathematical model that incorporates cell-gel interactions together with osmotic pressure to study the mechanical behaviour of biological gels. In particular, we consider an experiment where cells are seeded within a gel, which gradually compacts due to forces exerted on it by the cells. Adopting a one-dimensional Cartesian geometry for simplicity, we use a combination of analytical techniques and numerical simulations to investigate how cell traction forces interact with osmotic effects (which can lead to either gel swelling or contraction depending on the gel’s composition). Our results show that a number of qualitatively different behaviours are possible, depending on the composition of the gel (i.e. its chemical potentials) and the strength of the cell traction forces. A novel prediction of our model is that there are cases where the gel oscillates between swelling and contraction; to our knowledge, this behaviour has not been reported in experiments. We also consider how physical parameters like drag and viscosity affect the manner in which the gel evolves.


Introduction
Biological tissues are composed of cells living in extra-cellular matrix, hereafter designated ECM (Iordan et al. 2010). The ECM provides mechanical support to the cells in vivo and helps to regulate cell behaviour, as well as playing a key role in the mechanical behaviour of the tissues themselves (Humphrey et al. 2014;Dolega et al. 2021). The ECM in vivo can be thought of as a fibrous polymer network; it can consist of a number of different substances, including polysaccharides, proteoglycans and proteins. To reproduce in vivo-like behaviour when cells are cultured in vitro, they are often seeded in a gel, which aims to mimic the ECM. Since the structural protein collagen is the primary component of the ECM in many animal tissues, collagen gels are frequently used in laboratory studies, but a wide range of other natural (e.g. Matrigel) or synthetic (e.g. poly(lactic acid)) gels are also used (Wade and Burdick 2012). Improved understanding of the mechanical behaviour of biological gels, together with cell-cell and cell-gel interactions, will lead to better understanding of the development and functioning of tissues.
The mechanical characteristics of a tissue can have a powerful effect on cell behaviours such as proliferation, differentiation and cell motility; an effect that is, in fact, reciprocal, since the tissue is maintained by these cells (Rozario and DeSimone 2010;Wade and Burdick 2012;Humphrey et al. 2014). Experiments aimed at gaining more insight into the cell-ECM relationship and how each regulates and affects the other are often conducted using cell-seeded gels grown in vitro, where it is often easier to control conditions and make observations. One such experiment, presented by Moon and Tranquillo (1993), studies the contraction of collagen gels under the influence of cell traction forces. The cells interact mechanically with the polymer network surrounding them, leading to this polymer network being reorganised and compacted. These processes of ECM remodelling are important in tissue growth and development and, accordingly, are important in a range of related topics such as wound healing and tissue engineering (see e.g. Stevenson et al. (2010); Dyson et al. (2016) and references therein).
In the Moon and Tranquillo (1993) experiment, a microsphere of collagen gel is prepared and cells are seeded within the gel. The gel is placed in a bath of solvent (which provides nutrients for the cells) and these cells, through the traction stresses they exert on the polymer network, compact the gel over approximately two days. Although the degree of gel contraction can provide a measure of the forces exerted by the cells, it is dependent on the specific procedures employed in their experiment (e.g. cell density, gel composition and size, etc.). To obtain a quantitative measure of the cell-derived forces, Moon and Tranquillo (1993) developed a mathematical model of the process based on the mechanochemical theory of Murray et al. (1983).
The Moon and Tranquillo model assumed that the only forces acting on the gel were those exerted by the cells; however, gels can swell or contract in the absence of cells, for example, due to osmotic effects whereby solvent molecules can enter or leave the gel as a result of differences in the chemical potentials across the gel-solvent interface (Hong et al. 2010). This is often studied mathematically using multiphase flow models, in which the behaviour of both the polymer and solvent components of the gel are each accounted for. Recent studies have suggested that these osmotic effects could be used to manipulate the mechanical environment of cells in vitro by, for example, applying a compressive force to the gel in which they reside (Monnier et al. 2016).
Our aim in this paper is to gain more understanding of how cell and osmotic forces interact within cell-seeded gels grown in vitro through the development and investigation of a mathematical model. In particular, we are interested in the different qualitative outcomes that may arise (e.g. gel contraction, swelling or dissolution), and how these outcomes, together with the dynamics of the process, depend on factors such as the gel composition, cell seeding density, cell traction strength and interphase drag between the network and solution phases of the gel.
The mathematical models of both Moon and Tranquillo (1993) and Green et al. (2013) treat the gel as a single material (either fluid or solid) which must be assumed to be compressible for it to be able to contract. However, many types of biological gel are largely made up of water, an incompressible fluid, bringing the appropriateness of this assumption into question. Therefore, in this paper, we investigate the possibility that gel contraction is a result of syneresis, that is, as a result of water being squeezed out of the gel. This possibility is noted in both Moon and Tranquillo (1993) and Green et al. (2013), but was not pursued further. Osmotically-driven movement of solution into and out of gels has been studied in the context of gel mechanics using Flory-Huggins theory (Kumar and Gupta 2003). This theory is derived from statistical thermodynamics, assuming that the polymer network and solvent which make up the gel exist on a two-dimensional lattice, and considering the number of ways in which the polymer and solvent molecules can be arranged on this lattice. From this, the entropy and enthalpy of the mixture is derived. Using these, we can calculate the free energy of the mixture per unit volume, and thereby, the osmotic pressure (Kumar and Gupta 2003;Winstanley et al. 2011). This term refers to the additional pressure that must be applied to reach an equilibrium between the gel and solvent separated by a semi-permeable membrane (Winstanley et al. 2011); it is thus a function of the difference in the chemical potentials between the gel and the surrounding solvent. Osmotic pressure will be included in our model through these chemical potential functions which are discussed further in Sect. 2.3. A number of previous studies have modelled the movement of solution into and out of gels and gel-like substances, usually in the absence of cells, for example biofilms (Cogan and Keener 2004;Winstanley et al. 2011), swelling polymer gels (Keener et al. 2011b, a), and general polymer solutions (Doi 2011;Doi and Onuki 1992). In the majority of these studies, the gel is modelled as a mixture of two interacting components: a polymer network, and a fluid solvent.
The multiphase or mixture theory approach provides an obvious framework for modelling the dynamics of gel swelling and contraction, and models of this type have been developed by Keener et al. (2011a, b) and Mori et al. (2013). They consider a gel consisting of two phases, polymer network and solvent, and present mass balance and momentum balance equations for each phase. Osmotic effects are incorporated in the momentum balance equations for each phase through terms involving spatial gradients of the chemical potentials. The chemical potentials depend on the Flory-Huggins free energy, and are functions of the volume fractions of network and solvent. We aim to combine their approaches with mechanochemical theory (Murray et al. 1983) to produce a new mathematical model for the mechanics of cell-seeded gels which incorporates osmotic effects.
The paper is structured as follows. In Sect. 2, we develop our multiphase model for gel-solvent interaction. This model includes the novel addition of cell-gel mechanical interactions, considering the cell traction stresses acting on the polymer network, alongside the osmotic effects introduced by free energy. As done by Keener et al. (2011a, b) and Mori et al. (2013) we non-dimensionalise the model and transform it into 1D Cartesian coordinates to obtain a simple model for analysis. In Sect. 4, the equilibrium conditions for the model are presented, as well as a short time solution. This short time solution allows us to evaluate how the model evolves away from steady-states and initial conditions and, so, to infer the stability of steady states. We implement the 1D model numerically in Sect. 5 to study the conditions under which swelling or contraction occur. Finally, in Sect. 6 we discuss the key results from this model and possible extensions for future investigation.

Mathematical model
We consider a gel seeded with cells, which sits within a surrounding bath of liquid solvent (e.g. nutrient medium). This gel-solvent system is sketched in Fig. 1. The domain is divided into two regions: the gel region g , and the surrounding region of pure solvent s . We note that the problem can be studied in different geometries (i.e. g need not be spherical). We let x denote position in and t denote time. The centroid of the gel is at x = 0 and the gel-solvent interface, denoted g (x, t) = 0, is the boundary between g and s . This interface between the gel and surrounding solvent can move over time with the movement of solvent between the two regions.
We use a multiphase modelling approach based on the work of Keener et al. (2011b) and Mori et al. (2013), modified to incorporate cells. The gel is assumed to be made up of two phases, polymer and solvent, each of constant density, with volume fractions denoted by θ p (x, t) and θ s (x, t) respectively. Hence, we define g to be the region where θ p > 0 and θ s > 0, and s to be that where θ p = 0 and θ s = 1. Cells are only present in the gel region g , since we assume they adhere to the polymer network of the gel. We note that in Iordan et al. (2010), gels were constructed with collagen concentrations ranging from 0.42 to 1.8 mg/mL and cell concentrations of = g ∪ s . g contains the cell population as well as positive volume fractions for both the polymer network and solvent, whereas s contains only solvent. g (x, t) = 0 is the moving boundary of the gel, also referred to as the gel-solvent interface. s is the external boundary of the domain 0.7×10 7 cells/mL, 1.17×10 7 cells/mL and 1.8×10 7 cells/mL. This provided cell volume concentrations respectively of 4.7%, 7.8% and 12%. In the Moon & Tranquillo experiments, gels were constructed with 3 mg/mL collagen concentration and seeded with much smaller cell populations of 5×10 4 cells/mL, 1×10 5 cells/mL and 3×10 5 cells/mL, indicating that the volume fraction of cells in these experiments would be even less significant. Hence, for simplicity, we assume that the volume the cells occupy within the gel is negligible; we therefore do not include a cell volume fraction and instead consider cell density n(x, t), where n(x, t) = 0 in s (similar to Barocas and Tranquillo (1994)). Thus, the no-voids condition is satisfied everywhere in the domain = g ∪ s . Moreover, the model presented below, while written for the gel region g , is also applicable to the solvent region s on setting θ p = n = 0 and θ s = 1.

Mass and momentum conservation equations
We assume the mass of both polymer and solvent is conserved (i.e. production and degradation of both species is neglected), so that where v p (x, t) and v s (x, t) are the polymer and solvent velocities respectively. Given the no-voids condition (1), adding Eqs. (2a) and (2b) yields where v = θ p v p + θ s v s is the volume-averaged velocity of the polymer-solvent mixture (Keener et al. 2011b); we choose to replace (2b) with (3). Note that in the solvent region s we simply have ∇ · v s = 0. For simplicity, cell proliferation and death are neglected, so the cell population is fixed. Cells are assumed to move by a combination of advection with the polymer network and unbiased random motion. In common with earlier models (e.g. Moon and Tranquillo (1993); Murray et al. (1983)) the latter effect is modelled by diffusion, with constant diffusion coefficient, D. Conservation of cells then gives ∂n ∂t We obtain equations for v p (x, t) and v s (x, t) by considering the momentum balance throughout the domain. This requires us to make appropriate assumptions about the mechanical properties of the cell and solution. Based on experimental data, Barocas et al. (1995) treated the gel as a Maxwell viscoelastic fluid. However, given the timescale of gel compaction, Green et al. (2013) noted that the Deborah number (which gives the ratio of elastic to viscous effects) found experimentally for gels like collagen is small (O(10 −1 ) − O(10 −2 )), meaning that elastic effects can be ignored. Hence, following Keener et al. (2011b) and Mori et al. (2013), we model both the polymer and the solvent phases as viscous fluids with a common pressure, P. The viscous stresses in the two phases are encapsulated by the deviatoric tensors σ p and σ s , where the polymer stress tensor σ p and the rate of strain tensor e p are defined by with the solvent tensors σ s and e s similarly defined (with subscript s in place of p). The constants η j and κ j ( j = p, s) are the dynamic and bulk viscosities of each phase i respectively, and I is the identity tensor. As in Keener et al. (2011b), we assume that the forces exerted on the two phases come from inter-phase drag (which is proportional to the product of the volume fractions of the two phases) and chemical potential gradients. In addition, we include traction stresses exerted by cells on the polymer network. Inertia can be neglected on the time and length scales typical of experiments such as Moon and Tranquillo (1993), so that the momentum balances for the two phases are given by In Eqs. (5a) and (5b), μ p (θ p ) and μ s (θ p ) are the chemical potentials for the polymer and solvent respectively, while ξ ≥ 0 is the constant drag coefficient. The traction force exerted by the cells on the polymer network is a novel addition in this context of multiphase gel modelling; it is incorporated as a body force acting on the gel and is given by the gradient of θ p G(n), where G is a scalar potential energy function (Mori et al. 2013). The forms of the cell potential function G and the chemical potentials μ p and μ s are detailed in Sects. 2.2 and 2.3 below.

Cell potential energy function
We assume the energy potential associated with the cells to be given by the Hill equation This differs from the function described in previous works (Murray 2001;Moon and Tranquillo 1993;Green et al. 2013) in having n 2 rather than n in the numerator. This means that ∂G/∂n > 0 for all n, which ensures that the cell traction force, acts in the direction of increasing cell concentration. The positive parameter τ 0 provides a measure for the strength of cell traction forces, and λ is a positive contact inhibition parameter, which reflects that the force exerted by cells decreases as the cells become more densely populated.

Chemical potentials
The chemical potential functions μ p and μ s describe the work done by the free energy in the polymer and solvent, respectively, to affect the swelling or compaction of the gel. These are defined as where f (θ p ) is the free energy per unit volume of gel (Keener et al. 2011b). The free energy function, derived from polymer physics, is defined below. The polymer chemical potential μ p describes the change in free energy resulting from an additional polymer unit being added to the gel, while the solvent chemical potential μ s describes the change in free energy from an additional solvent unit being added (Keener et al. 2011b). The free energy of the system per unit volume of gel is where k B is the Boltzmann constant, T is temperature, ν m is the characteristic volume of a monomer in our system, N is the chain length of the polymer, χ is the Flory interaction parameter and the constants μ 0 p and μ 0 s are dimensionless quantities known as the standard free energies of the polymer and solvent respectively. The logarithmic terms in the function describe the entropy of mixing polymer and solvent; these terms always encourage swelling in the gel. The latter terms involving χ , μ 0 p and μ 0 s can increase the tendency for the gel to swell or contract depending on the signs of these parameters. The χ term describes the energy of mixing, while the terms involving μ 0 p and μ 0 s describe the interaction energy in a pure polymer or solvent state respectively (Rubinstein et al. 2003).
In most of the relevant literature (e.g. Mori et al. (2013); Rubinstein et al. (2003); Zhang et al. (2008)) the standard free energy parameters μ 0 p and μ 0 s are not included, so that the free energy function represents only the interaction or mixing of the phases as opposed to the total free energy (Keener et al. 2011a). However, following the work of Keener et al. (see Keener et al. (2011b, a); Sircar et al. (2013)), we retain the standard free energy terms for generality. In the framework of Mori et al. (2013) that we adopt, we will see that these terms do not contribute explicitly to the final model, due to cancellations of terms involving f (θ p ) and its derivatives. They are, nevertheless, contained implicitly through the derivation of the mixing parameter χ (see Rubinstein et al. (2003) for further detail).
The Flory interaction parameter χ is a dimensionless parameter that characterises the nature of the interaction between the phases in the mixture: χ < 0 indicates attraction between the phases, and accordingly, mixing of these components being energetically advantageous; χ > 0 corresponds to repulsion between the polymer and solvent, resulting in the phases preferring to separate (Rubinstein et al. 2003).
As in Keener et al. (2011b), from Eqs. (8a) and (8b), we can derive further useful relations between the chemical potentials and free energy. Firstly, we have the relation, We also have that which indicates that at stationary points of f (θ p ), the chemical potentials μ p and μ s must be equal.

Initial and boundary conditions
To close our system of equations, we need to impose suitable initial and boundary conditions. The initial conditions for the volume fractions and cell density are given by where 0 < θ i (x) < 1 and n i (x) ≥ 0 for x ∈ p , and θ i (x) = n i (x) = 0 for x ∈ s . The initial gel-solvent interface is given by We take the centroid of the gel to be fixed in space and, therefore, have zero velocity at the origin for all time, while no slip and no penetration on the external boundary of the domain s is given by The gel-solvent interface g (x, t) = 0 moves over time due to movement of the polymer phase, so that its position is given by the kinematic condition We assume there is no diffusive flux of cells out of the gel at the interface, so that wheren is the unit outward normal vector on g = 0. The continuity of stress across g = 0 is described by where the prime denotes differentiation with respect to θ p (Mori et al. 2013 Finally, at the interface, we have where we have introduced the superscript e to clearly designate a quantity in the solvent domain s external to the gel. This condition describes how the difference in pressure, chemical potential, and solvent stress across the interface drives fluid flow into or out from the gel, at a rate proportional to the resistance R ≥ 0 of the boundary (see Eq. (3.15) in Mori et al. (2013)). We note that an increase in R increases the resistance of the boundary so that it is more impervious to solvent flow, while a decrease indicates that it is easier for fluid to move across the boundary in and out of the gel. With the resistance R = 0, the normal solvent stresses are equal to the pressure difference across the interface.

One-dimensional Cartesian model
For simplicity, in this paper we investigate the behaviour of a one-dimensional gel (denoting the Cartesian coordinate x) which is assumed to be symmetrical about x = 0. We also assume that all quantities are continuous and differentiable at x = 0. Thus, we consider 0 ≤ x ≤ L(t) as the gel domain g in which 0 < θ p (x, t) < 1, θ s (x, t) = 1 − θ p (x, t), and n(x, t) ≥ 0, while the polymer and solvent velocities are v p (x, t) and v s (x, t). There is a fixed symmetry boundary at x = 0 and a moving boundary at x = L(t) (equivalently, g (x, t) = x − L(t) = 0) on which the kinematic condition (16) becomes Outside the gel domain (x > L(t)), we have the solvent domain s with θ p e (x, t) = 0 and θ e s (x, t) = 1, where we have introduced the superscript e to clearly designate the solvent domain external to the gel. From hereon, this superscript notation will be used for all quantities in s , while lack of the superscript e denotes quantities in g .
Since, by symmetry, v p (0, t) = v s (0, t) = 0, the continuity condition (3) implies that throughout g we have Similarly, throughout s we simply have v e s = 0, which satisfies the no-penetration condition (15) on the boundary s .
In the gel domain g , from the mass conservation Eq. (2a) we also have while the advection-diffusion equation for cell density (4) becomes and the momentum Eqs. (5a) and (5b) are now On multiplying (25) by θ s and (26) by θ p and taking the difference, we eliminate the pressure terms from the momentum equations, yielding On substituting for v s and using Eq. (28), (27) becomes We use (29) to replace (25).
In the solvent region s , where v e s = 0, θ p e = 0 and θ e s = 1, the solvent viscous stress tensor σ e s is zero, and from (8b) and (9), where f (0) is a constant. From the definition of f (θ p ) in Eq. (9), we see that f (0) = μ 0 s . In s , the momentum Eq. (5b) therefore simplifies to and we see that P e is at most a function of time t. The 1D form of the interface condition (19) while (20) Using (33) to eliminate pressure from (32) and using (22) to substitute for v s , we obtain From (14), the velocity at the origin is zero for all time, Given that v p and ∂v p /∂ x are continuous and differentiable at x = 0 and that v p is an odd function, it can be shown from Eq. (23) that if ∂θ i (0)/∂ x = 0, then ∂θ p (0, t)/∂ x = 0 for all t, i.e. if the polymer fraction θ p is initially symmetric and continuous about the origin, then it will remain so for all time. Since we take θ i (x) to be differentiable with ∂θ i (0)/∂ x = 0, we therefore have Finally, the symmetry of the cell density at x = 0 requires The no-voids condition (1), the conservation Eqs. (22), (23), (24), the momentum Eq. (29), and boundary conditions (21), (34), (37), (38) comprise a complete model for the polymer and solvent volume fractions θ p (x, t), θ s (x, t), the cell density n(x, t), and the polymer and solvent velocities v p (x, t), v s (x, t) in the gel, along with the length L(t) of the gel. To solve for the pressure difference between the gel and solvent regions, P(x, t) − P e (t), we must add the momentum equation (26) and the boundary condition (33). However, we choose not to solve for the pressure in the gel, given that we can study the mechanics driving the gel without its inclusion, and so drop (26) and (33) from the model.

Non-dimensionalisation
Let L = L(0) be the length scale, N be a characteristic cell density (typically the mean initial density), and let the time scale be the ratio of polymer viscosity η p to the free energy scale k B T /ν m . Using these scales, we non-dimensionalise our model variables as follows, where tildes denote dimensionless quantities, We also define the following dimensionless model parameters, again denoted by tildes, and the dimensionless free energy, chemical potential functions, and cell potential energy, The mass balance Eqs. (22), (23) and (24) are unchanged in form on writing them in terms of the scaled variables and parameters. Similarly, (8a), (8b), (10), (11) and the boundary conditions (21), (37), (38) are unchanged in form. Hence we do not re-write them here. The scaled forms of the momentum equation (29) and interface stress condition (34) become over 0 ≤x ≤L(t) and, atx =L(t), We now introduce a change in coordinates to shift our moving boundary problem onto a fixed domain. We define new coordinates X =x/L(t) and T =t, so that the domain 0 ≤x ≤L(t) is mapped to the fixed domain 0 ≤ X ≤ 1. On this fixed domain the model becomes, using dots to denote differentiation with respect to time T , primes to denote differentiation with respect to θ p , and dropping tildes on dimensionless variables and parameters for convenience, with boundary conditions at and boundary conditions at X = 1, In addition, we must specify suitable initial conditions with θ i (X ) chosen such that ∂θ i (0)/∂ X = 0. This completes our derivation of the 1D model, which comprises Eqs. (41a)-(41e) for θ s , θ p , n, v s and v p , together with (43c) for the gel length, L.

Steady state and short-time behaviour
We now consider steady state (i.e. long time) and short time solutions of our model. The former allows us to understand the necessary conditions for the gel to equilibrate. The latter provides insight into the stability of steady states, and helps us to verify our numerical simulation methods.

Steady state conditions
The system reaches equilibrium when θ p and n are such that there is zero net force everywhere, the velocities of polymer and solvent are zero everywhere, and . L = 0, i.e. the free boundary stops moving. We find that, for θ p and n, both spatially uniform and non-uniform steady state solutions are possible.
At a steady state, the mass and momentum conservation Eqs. (41c)-(41e) are To demonstrate that the velocities are zero at equilibrium, we integrate the steady state polymer advection Eq. (45a) with respect to X and apply the boundary condition v p = 0 at X = 0. This gives θ p v p = 0, and since we must have θ p > 0, we see that v p = 0 must hold at equilibrium. Accordingly, we must also have v s = 0 using Eq. (41b).
The momentum balance Eq. (45c) then yields the following equilibrium condition within the gel, From (8b), we note that and accordingly, Eq. (46) can be expressed in the form i.e. for the gel to be in equilibrium, the cell traction force must be balanced by the force due to chemical potential gradients. We note that with n = 0 (and hence G = 0), this is the same condition as in Keener et al. (2011b). Equation (47) is subject to the condition (43a) at the interface X = 1, which, at equilibrium, gives After integrating (47) and using (48) to set the constant of integration, we obtain which applies everywhere at equilibrium. We note that, as shown in (30), μ e s = f (0) is constant.
From the cell advection-diffusion Eq. (45b), we see that, at equilibrium, we must have For D = 0, this is trivially satisfied, and Eq. (49) is sufficient for the gel to equilibrate. In this case, it is possible to have equilibrium solutions where θ p and n depend on X .
With D = 0, after integrating (50) and applying the no-flux condition (43b), we find that i.e. at equilibrium, n must be spatially uniform. Now, given spatially uniform n, Eq. (46) can be written indicating that either G − θ p f or ∂θ p /∂ X must equal zero for equilibrium. Given the functional form of f (θ p ) (see (9)), we have Evaluating G − θ p f = 0, we find the following quadratic expression in θ p , This shows that, at most, we can have two positive solutions for θ p , depending on the values of the model parameters. However, given that θ p must be continuous, only a constant value of θ p will satisfy this condition. Thus, we must have ∂θ p /∂ X = 0 to satisfy G − θ p f = 0 at equilibrium. Therefore, from Eq. (52), we see that we must have at equilibrium, i.e. spatially uniform θ p . Therefore, if diffusion D = 0, n and θ p must be spatially uniform and satisfy Eq. (49) for the gel to reach a steady state. We note that, since the total masses of polymer and cells are conserved, any changes to the initial polymer fraction or cell density that change the total polymer or cell mass will lead to changes in L, θ p and n at later points in time. Hence, the equilibrium solutions for θ p , n and L are not independent of the initial conditions. We can use Eq. (49) to calculate the spatially uniform equilibrium values of θ p and n for a given set of parameter values (note that this does not preclude the existence of spatially-varying steady states for such parameter values when D = 0). This is useful for two reasons: firstly, it provides a method to confirm that numerical simulations (such as those we will see in Sect. 5) find the correct equilibrium values; secondly, it allows us to analyse how the equilibrium values of polymer and cell density change as chosen parameter values are adjusted. Together with a condition derived from the short time solutions in Sect. 4.2.2, we will use this in Sect. 4.2.3 to analyse steady state values of particular variables and parameters, as well as the stability of these equilibria as different parameter values change.

Short-time analysis
We next study the behaviour of the system (41a)-(44) on a short time scale to investigate the early time evolution of the system from non-equilibrium initial conditions. In Sect. 4.2.2, we will determine how the system evolves over small time in response to small spatial perturbations to equilibria.

Evolution from non-equilibrium initial conditions
We proceed by introducing the short time scale δ 1 and define T = δT . We then write the dependent variables as power series in δ, expanding about the initial conditions: with expansions for θ p and n similar to that for v p . Here L 0 = 1, v 0 (X ), θ 0 (X ), and n 0 (X ) are the initial conditions. For simplicity, and in the interests of finding an analytic solution, we shall restrict our attention to spatially uniform initial conditions for θ p and n, i.e. θ 0 (X ) = θ i and n 0 (X ) = 1, where θ i is constant and we have scaled n using the initial cell density as its characteristic value.
We substitute these expansions into (41c)-(43c), and solve to obtain: where we have introduced the parameters From these solutions we see that the evolution of the gel from uniform initial conditions depends on the sign of A 0 , which is determined by the first term in brackets in (59). Hence, the expansion or contraction of the gel depends on the balance between the initial cell and chemical potentials. With positive A 0 (e.g. with a large negative mixing parameter χ in μ s encouraging gel swelling), the length L increases as the gel stretches in response to the influx of solvent, causing θ p and n to decrease. Conversely, for negative A 0 (e.g. driven by a large cell traction parameter τ 0 ), the gel contracts in length as solvent is forced out, with θ p and n increasing accordingly. (Note that (49) implies the relevant term, and accordingly A 0 , is equal to zero if the initial condition θ 0 = θ i , n 0 = 1 is a steady state.) Since ξ > 0, then α > 0 and the cosh(α X ) functions in (57b) and (57c) increase monotonically with X so that θ p and n change most rapidly at the gel's interface. We note that, in the limit where there is no drag (ξ = 0), the small-time solutions corresponding to (57b) and (57c) are spatially uniform.
These solutions show us how the early time expansion or contraction behaviour of the gel depends on both the chemical and cell potentials when the gel is not initially in equilibrium. In the following section, we will consider how the gel behaves in response to spatial perturbations of a steady state, with the aim of gaining insight into the stability of equilibria.

Short-time behaviour of spatial perturbations to equilibria
We now examine how the system evolves over short time from initial conditions that are small amplitude spatial perturbations to equilibrium solutions. This will suggest the stability of the equilibrium state: an equilibrium will be taken as unstable if spatial perturbations increase in amplitude over time, leading the system to evolve away from the equilibrium; an equilibrium will be taken as stable if the perturbations decay. (We note that these stability criteria are supported by our numerical solutions to the model.) For simplicity, we restrict our attention to spatially uniform equilibria as required in the general case where D = 0.
We denote the dimensionless steady state by asterisks, L * , θ * , n * , v * , where v * = 0, and L * = n * = 1 (since length and cell density are scaled on their equilibrium values). We take δ to be the short time scale as in the previous section and let be the amplitude of the spatial perturbation, where δ 1. Next, we take the series (56), etc., and expand each of the terms L j , v j , θ j , n j , j = 1, 2, . . ., in powers of with the initial conditions We set θ 01 = cos(γ X ), n 01 = N 01 cos(γ X ), where N 01 is an O(1) constant describing the magnitude of the spatial perturbation to the cell density. We note that higher order terms of v 0 must be determined in such a way as is consistent with the other initial conditions. Note that θ 0 and n 0 satisfy the symmetry boundary conditions (42) at X = 0 for any choice of γ , while the no-flux cell boundary condition (43b) at X = 1 requires that γ = Z π for some integer Z . For Z = 0, the spatial perturbation is constant and so effectively only shifts our initial condition, resulting in a similar solution to that presented in Sect. 4.2.1. Furthermore, changing the sign of Z does not change θ 0 or n 0 . We therefore restrict our analysis to positive values of Z . We also note that the condition γ = Z π ensures that our choices of θ 0 and n 0 are such that the total masses of polymer and cells over the domain 0 ≤ X ≤ 1 are unchanged from the unperturbed initial masses (θ * and 1, respectively). Incorporating these spatial perturbations in , the expansions (56), etc., become and so on. We substitute these into the governing equations to obtain solutions describing the small time behaviour. Since the calculations are standard but somewhat lengthy, the details are omitted (but can be found in Reoch (2020)). The solutions for the case with non-zero drag (ξ = 0) are: where α is as defined in (58) on replacing θ i with θ * , and we have introduced the parameters For the no drag case (ξ = 0) the solutions simplify to: We note that, at O(δ ), Eq. (60b) does not satisfy the no-flux cell boundary condition at X = 1 due to the cosh(α X ) term. As explained in (Reoch 2020), this is because we neglect a higher-order term involving ∂ 2 n 11 /∂ X 2 , meaning that we have a singular perturbation problem and cannot satisfy all boundary conditions for n 11 . For the purposes of our analysis herein, we will continue to discuss this solution, as the error will be confined to the small region near X = 1. We also observe that in the zero-drag case, the boundary condition at X = 1 is, in fact, satisfied.
We see from Eq. (60c) that L increases or decreases depending on the sign of A 01 , which in turn depends on the value of z and the choice of γ . The A 01 sinh(α) term describes the 1D expansion or contraction of the gel over the short time scale as a result of the spatial perturbation. The changes in θ p and n due to this expansion or contraction are described by the A 01 cosh(α X ) terms in their solutions; as previously observed, cosh(α X ) increases monotonically with X since α > 0. Greater changes in θ p and n therefore occur as X increases across the spatial domain. This is similar to the solution given by Eqs. (57a)-(57c), where the gel length is governed by a sinh(α) term, while the cosh(α X ) terms determine the increase or decrease of the polymer fraction and cell density.
The trigonometric terms in the solutions for θ p and n, (60a) and (60b) respectively, describe whether the initial spatial perturbations increase in amplitude or decay over time. Growth in these perturbations implies an unstable equilibrium where the gel evolves away from its steady state; decaying perturbations meanwhile imply a stable equilibrium. In the zero-drag solution (64a)-(64c), L remains constant to O(δ ) and no hyperbolic functions appear. Accordingly, θ p and n evolve in space with no change in the gel's length, i.e. there is no flow at X = 1 when ξ = 0. Thus, a change in gel length only occurs over short time when ξ > 0. Since we are primarily interested here in the evolution of the polymer and cell distributions, we focus our attention on the case where ξ = 0. In this case, changes in the amplitude of the perturbations are simply governed by the coefficients of the trigonometric terms, avoiding any complications resulting from the small changes in gel length with drag present.

Steady state stability conditions
We now use our results from Sect. 4.2.2 above to investigate the response of steady state solutions to spatial perturbations. As explained above, to simplify the interpretation of the stability results presented in this section, we restrict our attention to the case where ξ = 0. By considering the amplitudes of the cos (γ X ) terms in (64), we see that for the amplitude of the initial perturbations to both θ p and n to be decreasing in time, and accordingly, for the solution to revert back to equilibrium, we require which can be expressed in the form Conversely, in the case that z < 0, the amplitude of the perturbation to θ p will grow with time. Hence, the growth or decay of perturbations in this system is dependent on the balance between free energy and cell force. In the absence of cells, it is the sign of f (θ * ) which determines the stability of an equilibrium. Given that where N is typically large, we see that it is the sign and magnitude of the mixing parameter χ together with the equilibrium fraction of solvent 1 − θ * that primarily determine this. Adding cells to the model will always reduce the value of z, and with other values held constant, move the equilibrium towards an unstable state. Similarly, with cells present, increasing cell traction strength τ 0 or reducing contact inhibition λ will make instability more likely.
Using the stability condition (66), we can infer whether equilibria are stable or unstable and study how both the equilibrium and its stability change as we adjust particular parameter values. The diagrams presented in this section have been generated by solving the equilibrium condition (49), which, upon substituting for G and μ s , and recalling that we have scaled the system such that the equilibrium cell density n * = 1, becomes while the first term vanishes for n * = 0. We solve this expression for a chosen parameter as we vary θ * between 0 and 1, holding other parameters fixed. Alternatively, we find the relationship between two parameters using (68) while keeping θ * and other parameters fixed. Through this analysis, we can also determine whether spatially uniform steady states exist for a particular set of parameter values and initial conditions. Figure 2a-c demonstrate how θ * varies as we change the mixing parameter χ at different values of τ 0 (we note that changes in the dimensionless parameter τ 0 can correspond to changes in the characteristic cell density-here the physical steady state value-or the dimensional cell traction strength). In the majority of cases, larger values of χ indicate greater levels of contraction in the gel, corresponding to larger values of θ * ; this outcome should be expected, as increasing χ indicates that separation of the two phases in the gel is more favourable. In Fig. 2a, it is shown that in the absence of cells (n * = 0), two equilibrium values of θ * -one stable and one unstable-exist for the same parameter values. Note that the stability of these steady states is inferred using Eq. (66). We also see that, in the absence of cells, steady states θ * exist only for positive values of χ (given N = 100); with χ < 0, the terms in the free energy function all promote mixing between solvent and polymer, and accordingly, the gel keeps expanding until it dissolves. In this example, if a gel's initial fraction of polymer θ i is in the region beneath the blue solid line, the gel will contract to equilibrium with a greater value of θ p , while if θ i is above the branch of stable equilibria, the gel will swell to a steady state. For χ < 0.62, there are no steady states possible for any initial condition θ i (i.e. the gel dissolves). This indicates that small changes to the initial composition of a gel, e.g. the fraction of polymer or make-up of the solvent, could have significant impacts on its subsequent behaviour and possible steady state.
With cells introduced into the system, θ * increases monotonically as χ increases, as seen in Fig. 2b and c where τ 0 = 0.25 and τ 0 = 1 respectively. In these examples, we see that there are no longer any unstable equilibria, and that stable equilibria now exist over the spectrum of χ values. We see that as the traction parameter increases from τ 0 = 0.25 (Fig. 2b) to τ 0 = 1 (Fig. 2c), the equilibrium polymer fraction is greater for the same values of χ . For example, in Fig. 2b where τ 0 = 0.25, at χ = 0 we have θ * = 0.2, while in Fig. 2c where τ 0 = 1, at χ = 0 we have θ * = 0.58. Similarly, we see in these figures that with τ 0 increasing from τ 0 = 0.25 to τ 0 = 1, the same particular equilibrium value θ * is found with a decreasing value of the mixing parameter χ .
This relationship between τ 0 and χ is reinforced in Fig. 2d, where χ is plotted against τ 0 for fixed θ * = 0.5. We see that χ decreases linearly with increasing τ 0 ; as the cells exert more force, lower values of the interaction parameter are needed to keep the system at the same equilibrium value of polymer. Similarly, with a larger value of χ , less cell traction is necessary to maintain this equilibrium. This linear relationship between χ and τ 0 when the polymer fraction is fixed can be clearly seen in Eq. (68).
In Fig. 2e it is shown that, as would be expected, for fixed χ , larger values of τ 0 correspond to larger θ * , i.e. greater compaction in the polymer network. We note that we assume τ 0 ≥ 0; therefore, the equilibria that cross the vertical line at τ 0 = 0 are not relevant to the current study. However, plotting for τ 0 < 0 reveals a very small branch of permissible unstable equilibria in the region approaching θ * = 0, but the vast majority of initial conditions here will reach a stable steady state. In the small region where two steady states exist, gels with an initial polymer fraction above the unstable values of θ * will contract to the stable equilibrium, whereas those below the unstable values will swell until the gel dissolves.

Numerical simulations
We now perform numerical simulations to investigate the behaviours predicted by our model. We consider a range of initial conditions-both uniform and non-uniformand parameter values to better understand the emergent behaviours that can arise as a result of the interacting factors in the system. The one-dimensional governing Eqs. (41a)-(44) are simulated in MATLAB, using finite difference methods. Central differencing is used in the velocity Eq. (41e), excluding at X = 1, where a one-sided difference is used for derivatives of θ p , and the boundary condition (43a) is used to provide a ghost point for v p . The Crank-Nicolson method is used for Eqs. (41c) and (41d), with one-sided differences used for derivatives of v p at X = 1. We verified the code by comparing the simulation results with the short-time and steady state solutions  derived in § 4, and by checking the scheme conserves mass effectively over time. We calculated the percentage change in the mass of polymer and cells at the initial and final points in time for each of the simulations presented in this section: with a time step of dT = 0.0005 and spatial step of d X = 0.002, the worst-case change in mass for θ p or n was 0.0076%. Full details of the numerical scheme and verification checks can be found in Reoch (2020). The aim of our simulations is to illustrate the qualitative behaviours of our model for different initial conditions and in different parameter regimes. For all the simulations presented in this section, certain initial conditions and parameters are kept fixed: these are shown in Table 1. For example, in nondimensionalising the model, the length scale is set such that the initial length L(0) = 1. Similarly, when cells are present, we choose the average initial cell density as the characteristic value, so that n i = 1 for an initially uniform cell distribution. Given that μ 0 p and μ 0 s do not appear in the final set of model equations due to cancellation of terms when taking μ s − μ e s in interface condition (43a), we set these to zero. We set the bulk viscosities κ p and κ s to zero without loss of generality, as these terms only appear in linear combination with the dynamic viscosity parameters. The polymer chain length N is generally large for polymer and solution mixtures, therefore we set N = 100 (Rubinstein et al. 2003), and we set the contact inhibition parameter λ = 1. The parameters which are varied between simulations are listed in Table 2. Unfortunately, we lack experimental data which would allow us to obtain good estimates of these model parameters, and so in most cases we vary them over a range centred on one. We note that the mixing parameter χ is the only term appearing in the final system of equations for which negative values are relevant to this study.

Cell-free gel, uniform initial conditions
In a cell-free gel (where n = 0), swelling or contraction is driven by the free energy of the system; gradients in chemical potentials on either side of the gel-solvent interface induce the movement of solvent and polymer. This is similar to the results presented in Keener et al. (2011b). In Fig. 3a, where θ i = 0.6 and χ = 0.75, the balance in chemical potentials μ p and μ s produces an osmotic pressure gradient, causing solvent to enter the gel from the surrounding solvent region s ; the gel thus swells until an equilibrium is reached with θ * = 0.45 and L * = 1.34. Conversely, in Fig. 3b and c, we see the gel contract to an equilibrium state. The free energy in the system has been altered in two different ways here to induce contraction. In Fig. 3b, we have taken the same initial conditions as Fig. 3a, but with an increased strength of mixing parameter χ = 1.5. In Fig. 3c, the initial fraction of polymer has been decreased to θ i = 0.25, with the value of χ remaining at χ = 0.75. The effect in both instances is to increase the initial free energy in the gel, resulting in a situation where the gradient in chemical potentials will induce solvent to flow out from the gel to balance the potentials, and hence result in a smaller equilibrium length. The gel equilibrates with θ * = 0.86 and L * = 0.7 in Fig. 3b, and θ * = 0.45 and L * = 0.56 in Fig. 3c. These equilibria all clearly satisfy the mass conservation relation L * θ * = θ i (as indeed will all steady states found). We note that the simulations in Fig. 3a and c reach the same equilibrium value, θ * = 0.45, for the two different initial conditions; this corresponds to the equilibrium predicted in Fig. 2a with χ = 0.75 and the same fixed set of parameter values otherwise. Fig. 3b meanwhile confirms that, for the same initial conditions and parameter set, increasing the value of χ will result in an equilibrium with a larger polymer fraction (θ * = 0.86). This is also in agreement with Fig. 2a.
We also note that the polymer fraction at X = 1 (shown by red dashed lines) evolves slightly faster than that at X = 0 (blue solid lines). This lag reflects the time taken for the solvent to flow into or away from the centre of the gel. We discuss the parameters affecting this lag and the spatial profiles of the polymer as the gel evolves in Sect. 5.3.
We have shown here that, in the absence of cells, the gel will swell or contract depending on the balance between chemical potential gradients across the gel-solvent interface. These behaviours echo those found by Keener et al. (2011b), which is expected since our gel model builds on their work. Comparing our model to that of Keener et al., we note that while the inclusion of boundary resistance in our model only affects the rate of the gel's elongation or shrinking, the absence of any contribution from the standard free energy parameters μ 0 p and μ 0 s will change the final gel length and polymer fraction found here for the same set of parameters used in Keener et al. We next introduce cells into the simulations to study their effect on the gel's behaviour.

Cell-gel system
In Fig. 4a, we use the same gel parameters as for Fig. 3a, and introduce a cell population with weak traction (n i = 1, τ 0 = 0.1; note also D = 0.01). We see that the gel still swells to a steady state with the cell traction parameter set at this low level. However, compared to the simulation in Fig. 3a, the final size of the gel is now smaller (equilibrating here at θ * = 0.54, n * = 0.91, L * = 1.1), indicating that the cells are exerting some contractile force that counters the expansion due to osmotic effects. We increase the traction parameter to τ 0 = 1 in Fig. 4b; once cell traction is increased over a certain threshold, the gel will switch from expansion to contraction. In this instance, the cell traction stresses are stronger than the chemical potential gradient, and as the cells compact the polymer network, solvent is squeezed from the gel until it reaches a steady state once the mechanical forces are in balance, where θ * = 0.86, n * = 1.44, L * = 0.69. We have therefore established that introducing cells into a gel that would otherwise swell can induce a switch in behaviour, resulting in a significantly different outcome for the gel. As in Sect. 5.1, the equilibria found here satisfy the mass conservation relations L * θ * = θ i , L * n * = 1. The time taken to equilibrate is noticeably different across the simulations seen so far, as the rate at which the gel evolves is affected by the strength of a number of competing forces. For example, we see in Fig. 3a that equilibrium is reached at approximately T = 130, while in Fig. 3b, due to the larger value of χ , not only does the gel contract, but it equilibrates by T = 15. In a case like that presented in Fig.  3b, where the free energy alone induces gel contraction, adding cells to this gel will lead to a steady state being reached more quickly (result not shown). Therefore, the magnitude of parameters like the interaction energy χ and cell traction τ 0 will affect the time taken to reach a steady state. Alongside this, mechanical factors like drag and viscosity will impact the gel's temporal evolution.

Effects of mechanical parameters and diffusion on gel evolution
We now study how the ratios of drag ξ , resistance R, and solvent viscosity η s relative to polymer viscosity η p affect the rate at which a gel evolves to equilibrium and the manner in which it does so spatially. In the simulations presented in Figs. 5a, 6c, we take a gel with the same initial conditions, free energy parameters, and cell force parameters as that presented in Fig. 4b; this gel will therefore reach the same equilibrium regardless of parameters like drag and viscosity (given that the equilibrium is determined by Eq. (49), in which these mechanical parameters do not appear). We first set ξ , R and η s to be small, corresponding to the situation where these effects are insignificant relative to polymer dynamic viscosity, and as such, polymer dynamic viscosity is the dominant mechanical characteristic; we will refer to this as the base case in the comparisons that follow. We note that due to the scaling used here, we have effectively set η p = 1. To better understand the impact of each parameter on the gel's evolution, we then change one of ξ , R and η s in turn while holding the others constant. We take D = 0 here so that diffusion has no impact on the cell and polymer distributions. The gel reaches the same steady state in each case (θ * = 0.86, n * = 1.44, L * = 0.69), albeit at different times and with different lags between X = 1 and X = 0.
In Fig. 5a, we compare the evolution of θ p for the base case with ξ = R = η s = 0.1 (shown by the red solid line for X = 0 and the red dashed line for X = 1) and for a gel with large drag, where ξ = 5 and R = η s = 0.1 (shown by the blue solid line for X = 0 and blue dashed line for X = 1). We see that increasing the drag coefficient slows down the evolution of the polymer fraction (with the gel equilibrating at T ≈ 25 with large ξ compared to T ≈ 10 with small ξ ). Furthermore, for large drag, θ p changes at a much slower rate at X = 0 than at X = 1, while there is little difference in θ p between X = 0 and X = 1 when drag is small. This is reflected in Fig. 6a and b, which show the spatial profiles for θ p at increasing points in time for the base case and large drag case respectively. With polymer viscosity dominant in Fig. 6a, the gel evolves across the spatial domain in a largely uniform manner. In Fig.  6b, much stronger spatial variations are evident, reflecting that, while the gel evolves quickly at the interface, it takes much longer for solvent to flow through the domain due to the extra resistance when drag is large. In practice, the drag will depend on properties of both gel and solvent-e.g. the permeability of the polymer gel used, and the viscosity of the solution. It could be varied by changing the culture medium or the type of polymer used to make the gel. However, we note that experimentally, it would be difficult to change either of these parameters independent of others in the model (e.g. changing the type of solvent would change both drag and the solvent viscosity).
In Fig. 5b we similarly compare the gel's behaviour with interface resistance R and solvent viscosity η s each large relative to polymer viscosity. For large resistance (shown by the green solid and dotted curves for X = 0 and X = 1 respectively), we take R = 5, ξ = η s = 0.1, while for large viscosity (shown by the black solid and dotted curves for X = 0 and X = 1 respectively), we set η s = 5, ξ = R = 0.1. Note that this figure should be compared to the base case given by the red curves in Fig. 5a. Increasing the resistance parameter R, the rate of change of θ p is slowed at X = 1 compared to the base case (with equilibrium now reached at T ≈ 20). This reflects the fact that the boundary of the gel is less permeable to fluid flow with larger R. In this case, the polymer fraction remains almost uniform across the spatial domain, i.e. there is no additional lag induced between X = 1 and X = 0 with large R. Large solvent viscosity η s has the effect of slowing down gel contraction further still, with the gel not equilibrating until T ≈ 80 (note that equilibrium for the large η s case is not shown in Fig. 5b). As with R, in this case there is minimal lag between the evolution at X = 1 and X = 0.
We now study the effect of non-zero diffusion on the gel's behaviour. The diffusion coefficient models the degree of random motion of the cells, which could be varied experimentally by e.g. using different cell types with differing degrees of motility, or by adding chemicals that promote undirected cell movement (e.g. chemokines) to the culture medium. As the diffusion coefficient increases, we expect to see more uniform spatial profiles in the cell density as well as the polymer fraction, as cells spread more evenly across the gel through random motion. In Fig. 6b we showed the spatial profiles of θ p for a gel with a large drag coefficient and zero diffusion. In Fig. 6c, we see that the spatial distribution of cells in this gel is similarly non-uniform over much of the gel's evolution. We now take this same gel with large drag, but introduce cell diffusion, setting D = 0.005; the time evolution in this case is shown in Fig. 7a. In Fig. 7b, we see that with large drag, the cell density initially increases rapidly at X = 1 (like in Fig. 6c). Over time, cells move down their density gradient towards X = 0 due to the diffusion term now present. The cells become more dense in this region, eventually overshooting the equilibrium value (see n in Fig. 7a); however, as time progresses, the diffusive movement leads n to converge to its equilibrium value across the domain.
Increasing diffusion further to D = 1 in Fig. 7c, the strength of the random motion is such that the cells remain well spread spatially for all time. We note that in this case, the polymer still evolves with a lag across the spatial domain, like seen in Fig.  6b (result not shown).

Reduced initial polymer fraction
Taking a smaller initial polymer fraction, we can see different dynamics emerge in the evolution of the gel. Setting θ i = 0.2 and all other parameters as in Fig. 4b, we see in Fig. 8a and b that the polymer fraction and cell density evolve slowly over the beginning phases (T = 0 to T ≈ 2), before a period of rapid increase (T ≈ 2 to T ≈ 7), which slows down again as the gel moves towards its steady state (T ≈ 7 to T ≈ 12). From this lower initial value of θ p , much greater contraction is evident in the gel, which  Fig. 4b, we saw that with θ i = 0.6, the resulting steady state was (θ * , n * , L * ) = (0.86, 1.44, 0.69). This demonstrates that, with the initial cell density constant, there is a negative correlation between the initial fraction of polymer in the gel and the degree to which both the gel contracts (seen in L * ) and the polymer compacts (seen in θ * ). The reason is that for spatially uniform steady states and initial conditions, by mass conservation we must have L * θ * = θ i , L * n * = 1; this implies that n * = θ * /θ i . We see that the decrease in θ * and n * with θ i increasing, as in the examples above, requires an increase in the equilibrium length L * to conserve mass for these quantities. This negative correlation was also observed in experiments presented in Stevenson et al. (2010). This example is discussed further in relation to other models in Sect. 6. We note that we see similar behaviour take place in a cell-free gel. For example, amending the example in Fig. 8a such that θ i = 0.2, n i = 0, χ = 1.5, the gel evolves in a similar manner, with a slow initial phase of polymer compaction followed by rapid evolution and then slowing near the steady state (see Fig. 8c). In this instance, the gel equilibrates with θ * = 0.86, L * = 0.23. Furthermore, taking θ i = 0.4 and θ i = 0.6, we reach steady states of (θ * , L * ) = (0.86, 0.47) and (θ * , L * ) = (0.86, 0.70) respectively (result not shown for θ i = 0.4; see Fig. 3b for θ i = 0.6). While the equilibrium fraction of polymer is the same regardless of the initial condition in the cell-free case (provided the same free energy parameters are used), we see that the change in gel length is also negatively correlated with the initial polymer fraction here. This is clear from the mass conservation condition L * θ * = θ i , given θ * is fixed by the free energy parameters in the cell-free case. These results indicate that the presence of cells is therefore necessary to see the negative correlation between initial and final polymer fractions.

Non-uniform initial conditions
The numerical simulations presented thus far have been performed using spatially uniform initial conditions. Despite spatial variations being evident while the gel is evolving, these initial conditions eventually produce spatially uniform steady states. This matches previous work such as Keener et al. (2011b) looking at cell-free models, wherein only spatially uniform equilibria are found. We now consider examples with spatially non-uniform initial conditions, finding that these initial conditions can result in spatially varying equilibrium solutions. This is a novel behaviour arising in our model from the presence of cells. We will consider non-uniform initial conditions in both the polymer and cells separately.  2.5, 5, 10, 20, 30, 40, 150 We first evaluate a cell-free gel with a spatially varying initial polymer distribution; this allows us to establish a baseline against which we can evaluate the impact of cells. We take a gel as specified in Fig. 3a where θ i = 0.6 and χ = 0.75, for which the gel swells to an equilibrium with θ * = 0.45 and L * = 1.34. We add a spatially varying component to the initial condition for the polymer fraction, such that θ i = 0.6 + 0.025 cos(π X ); this corresponds to a gel where the polymer is slightly bunched at the gel's centre (where X = 0) and less than the mean value at the gel's edge (where X = 1). We note that this initial condition satisfies the symmetry condition ∂θ p /∂ X = 0 at X = 0. In Fig. 9a, we see that this gel swells to the same steady state as for uniform θ i = 0.6, and evolves on a similar time scale (see Fig. 3a for comparison). In Fig. 9b, we see the spatial distribution of polymer across the length of the gel at increasing points in time; the polymer, initially more concentrated towards X = 0, smooths out over time as the gel swells, eventually becoming uniform as it expands to its steady state where θ * is constant. Therefore, in the simulations we have seen, spatial variations in the initial polymer distribution in a cell-free gel do not affect the equilibrium outcome.
We now take the same gel with a spatially varying polymer initial condition and include a cell population where n i = 1 and τ 0 = 1, noting that D = 0. The time evolution for this gel is shown in Fig. 10a; we see that the gel, which swelled in the absence of cells due to osmotic effects, now contracts to an equilibrium with spatially non-uniform solutions for both the polymer fraction and cell density. The mean equilibrium values of θ p and n here, (θ * = 0.86,n * = 1.44), are the same as the steady state found in Fig. 4b. Figure 10b and c show the spatial distributions of θ p and n respectively over time. In contrast to the cell-free case, the cell forces present in the system are stronger than the chemical potentials, and so induce the gel to contract. We see that the polymer is initially less compacted towards X = 1; it must therefore contract more in this region to move to its steady state value. As the fraction of polymer increases in this region, cells then become more concentrated, which reinforces this The presence of non-zero drag also contributes to the formation of the spatial gradients seen, as discussed earlier. With D = 0, there is no requirement for the cells to even out over time, and therefore, we see a non-uniform distribution remaining at equilibrium. By the end of the process, the polymer fraction has ended slightly larger at X = 1 than at X = 0, reflecting the greater density of cells in that region. We next take the same system and change the spatial perturbation to the initial condition from the polymer fraction to the cell density, such that our initial conditions are n i = 1 + 0.05 cos(π X ), θ i = 0.6. This corresponds to a gel where the cells are now initially more densely seeded at X = 0 (i.e. the centre of a gel that is symmetric about the origin). In the time evolution for this system, shown in Fig. 11a, we see the gel reaches a steady state with the same mean polymer fraction and cell density as the previous example, albeit with different spatial distributions.  1,2,3,4,5,6,15) the polymer velocity v p goes to zero across the spatial domain over time, demonstrating that the gel reaches its equilibrium state. The spatial profiles here are displayed in Fig.  11c and d. The presence of drag creates a slight increase in both θ p and n at X = 1 while the gel evolves. As the gel approaches its steady state, the evolution at X = 1 slows while the polymer fraction and cell density continue to increase across the rest of the domain. Greater cell concentrations around X = 0 result in polymer being pulled to the gel centre, and finally, a higher polymer fraction in this region at equilibrium. At this resulting steady state, we see that the amplitude of the cell profile is greater than the amplitude of n i (the amplitude is approximately 0.04 at equilibrium and 0.025 initially), while the polymer fraction is slightly larger at X = 0 compared to X = 1, in contrast to the previous example.
These non-uniform equilibria have been found with diffusion D = 0. As shown in Sect. 4, for the gel to equilibrate with D = 0, n * and θ * must be spatially uniform. We have confirmed numerically that adding diffusion to these simulations will always result in a spatially uniform steady state on a long enough time horizon, with the additional random cell motion smoothing the cell profile, and subsequently, polymer profile as well (see (Reoch 2020) for examples).

Oscillating behaviour
Our model exhibits a novel behaviour where the cell density and polymer fraction will spatially oscillate as the system evolves, i.e. parts of the gel will switch back and forth between swelling and contraction over time; this is displayed in Fig. 12a. In this case, uniform initial conditions are used, with θ i = 0.5 and n i = 1, while there is a negative mixing parameter χ = −0.1 (indicating that the polymer and solvent would prefer to mix) and fairly strong cell traction τ 0 = 0.8. This behaviour emerges from the combination of parameters leading to the cell traction stress being finely balanced with the free energy. We have found that it also requires both the drag parameter ξ and the resistance parameter R to be sufficiently large (ξ = R = 1.5 here for example), otherwise these oscillations are not evident. We note that this behaviour occurs on a very long time scale and that the gel eventually dissolves in this situation (with θ p → 0). The length of the gel increases monotonically over time (as can be seen from Fig. 12b).
This behaviour comes about as a result of the competition between osmotic effects working to expand the gel and cell traction acting to contract it. The interface resistance slows the evolution at X = 1, while due to the presence of drag, steep gradients develop in the polymer fraction and cell density as the gel swells, with θ p and n decreasing most near X = 1. These significant gradients can be seen in the spatial profiles for θ p shown in Fig. 12b. Large variations in the cell density are similarly evident between the two ends of the spatial domain (result not shown). The greater density of cells around X = 0 produces a gradient such that the cell force starts to pull polymer back towards X = 0. The gel then contracts locally in this region, while still swelling across the domain towards X = 1. Eventually, the chemical potential gradients are such across the domain that the gel starts swelling for all X again, with these local fluctuations repeating once more as the gel slowly expands in a non-uniform manner. Over a long enough time frame, the gel eventually dissolves.
In Fig. 12c-e, we see the effect of reducing the interface resistance and drag in these simulations. With resistance and drag each taken separately to be small ( Fig. 12c and  d respectively), we see reduced oscillations in the gel; however, this behaviour still occurs. When both parameters are taken to be small, the oscillations are no longer present and the gel dissolves in a more typical manner (see Fig. 12e). Similarly, with very small diffusion (e.g. D ≈ 0.0001) the oscillations occur, but larger diffusion coefficients smooth out spatial gradients in the cell density and so prevent this oscillating behaviour from occurring (results not shown).
To our knowledge, this oscillating type of behaviour predicted by our model has not been reported experimentally. We note that it might not be easy to detect, since (as can be seen from Fig. 12b), the overall gel length increases monotonically in time. If we were to track a particular material point in the gel as it moves as the gel swells, we would notice an oscillation in the polymer volume fraction. We also note that the model only appears to produce such behaviour in a relatively restricted region of parameter space. However, if such behaviour were to be observed, it would help to validate our model.

Discussion
In this paper, we have presented a new model for cell-induced gel contraction, and studied its behaviour in a 1D Cartesian geometry. This has allowed us to develop a thorough understanding of the conditions under which the gel equilibrates, the conditions affecting the early time behaviour and the stability of the system, and, through numerical solution, the qualitative behaviours that can occur. Throughout, we have seen that the balance between chemical potentials and cell-derived forces is crucial to determining the gel's behaviour. We have shown that the presence of cells can cause a gel that would otherwise swell to contract; meanwhile, sufficiently strong osmotic forces can cause a gel to swell even with cells present. Moreover, the initial fraction of polymer was shown to negatively correlate with the final polymer fraction in cell-gel systems, and negatively correlate with the final gel length with or without cells present.
In Sect. 4, we studied the long and short time behaviour of the gel, showing that it is governed by the balance between cell and osmotic forces. The gel equilibrates when the cell and solvent potentials inside the gel are in balance with the external free energy; similarly, the early time behaviour and stability of equilibria are determined by the gradients of these functions inside the gel. Through deriving an analytical solution for the system's short time behaviour, we investigated how steady states respond (over short time) to spatial perturbations and determined whether these perturbations will grow or decay.
This analysis also allowed us to evaluate how steady state values and stability change with variations in parameter values. As discussed in Sect. 2.3, due to the boundary conditions we use, the standard free energy constants μ 0 s and μ 0 p do not appear in the final system of equations. As a result of this, we are not able to reproduce the examples of bistable equilibria presented in Keener et al. (2011a) where μ 0 p is varied against θ p . However, as mentioned in Sect. 2.3, the parameters μ 0 p and μ 0 s are typically not present in studies using the Flory-Huggins free energy.
In a laboratory setting, this analysis of parameter values and steady state outcomes can be used to predict experimental outcomes given specific gel configurations, e.g. suggesting whether a gel will equilibrate or dissolve, or if a different configuration is needed to reach a desired experimental result. This analysis may also allow for physical parameter values to be determined given comparison with experimental results. For example, if we are given experimental data for an initial gel configuration and its subsequent equilibrium state, we may be able to determine that particular parameters must lie within certain ranges through comparison with such steady state diagrams as presented in Fig. 2. Reducing the resistance parameter (to R = 0.5) or the drag (to ξ = 0.5) reduces the oscillation in θ p , as shown in c and d, respectively (note colour key and other parameter values are as for a). Reducing both (to ξ = R = 0.5) eliminates the oscillations altogether -see e In Sect. 5, we presented novel results relating to the gel's evolution, these being the presence of non-uniform equilibria and the case where the polymer fraction and cell density oscillate. Spatially non-uniform steady states were found to eventuate in the cell-gel system from non-uniform initial conditions in the polymer fraction or cell density in the absence of diffusion. In the oscillating case, due to competition between the free energy and cell traction, we see the fraction of polymer and the cell density repeatedly fluctuate within a spatial region as the gel swells, before it eventually dissolves. To our knowledge, this behaviour has not been described in the literature before and might be investigated experimentally.
Recent experimental work has suggested using osmotic pressure as a way to impose a desired mechanical compression on cells cultured in vitro (Monnier et al. 2016;Dolega et al. 2017). Our model provides a framework to quantify and evaluate such methods. Although to our knowledge, no one has yet used osmotic pressure to impose dynamic cycles of compression or tension on cells within a gel, our results suggest that this might be possible by, for example, changing the composition of the solution in the solvent bath surrounding the gel as a function of time. Again, this model could be used to predict the ensuing cycles of gel expansion and contraction, as well as to match the frequency and amplitude of these cycles to those seen in vivo. This might be beneficial in culturing cartilage cells for example, where oscillating stresses can lead to better mechanical and cell properties in the cells and tissues grown in vitro (Salinas et al. 2018). Oscillating fluid flow has also been seen to be an important mechanism in areas like proteoglycan production (Eifler et al. 2006) and regulating calcium concentrations (Edlich et al. 2001).
Our analysis has focussed on the qualitative behaviours predicted by our model in a simple 1D setting. To facilitate greater comparison with experiments, for example that presented in Moon and Tranquillo (1993), a couple of different steps could be taken. An extension to our work here, if a consistent set of data for relevant experiments was available, would be to fit such experimental data for our model parameters and initial conditions, allowing for a more direct comparison between these experiments and simulations like those presented in Sect. 5. Transforming the model to spherically symmetric coordinates would also help in comparing our results with models looking at spheres of gel like those in Moon and Tranquillo (1993) and Green et al. (2013), although we note that, given that such a model would be one-dimensional in the gel's radius r , we would not expect significantly different qualitative outcomes to those seen with our model here, which is one-dimensional in the gel's length.
In Fig. 8a and b, we saw the contraction of a gel with a small initial polymer fraction. In this instance, the polymer fraction and cell density increased gradually at the beginning and end of the gel's evolution, bookending a period of rapid contraction where θ p and n increased significantly. This behaviour is comparable to examples of mechanically-driven gel contraction presented in Moon and Tranquillo (1993) and Green et al. (2013). There was an initial lag in the evolution of the gel's radius seen in Green's model which replicates experimental observations from Moon and Tranquillo (1993). This initial lag was not present in Moon and Tranquillo's mathematical model, and similarly, we did not see an initial lag in changes to the gel's length in our simulations. Moon and Tranquillo posit that the initial delay seen experimentally is a consequence of the cells spreading after being seeded; it may not be present in our simulations as a consequence of the cells being smoothly distributed at initial time, therefore not requiring a lead time to redistribute themselves through the gel as may happen in vitro.
We also saw in Sect. 5.4 that, with all else held constant, the gel reached a smaller equilibrium length with a smaller initial polymer fraction. Similarly, with cells present, a smaller initial polymer fraction resulted in a larger value for the polymer fraction at equilibrium. Our model therefore captures the negative correlation between the initial polymer concentration and final concentration highlighted in the experimental study presented by Stevenson et al. (2010), who also reference this behaviour occurring in experiments such as Zhu et al. (2001) and Evans and Barocas (2009). In the absence of cells, in Fig. 8c and associated simulations, we saw that decreasing the initial polymer fraction θ i similarly led to gels with a smaller equilibrium length, i.e. gels that have contracted further. Without cells, the equilibrium value of the polymer fraction θ * is determined by the parameters in the free energy function; therefore, unlike what was seen for cell-gel systems, θ * remained the same with increases in θ i . We can therefore see that cell forces play an important role in the negative correlation between initial and final polymer concentrations seen experimentally.
Studies such as Barocas et al. (1995) that estimate the value of the cell traction parameter τ 0 typically use models that focus on cell-gel mechanics, i.e. they do not include the presence of chemical potentials. We have demonstrated herein that chemical potentials can counteract cell traction and affect the degree of gel contraction witnessed. This therefore indicates that models that neglect the free energy in the gelsolvent system may, in fact, underestimate the magnitude of cell traction stresses, since the degree of compaction in the experiment will also depend on the mixing energy of the polymer and solvent. It also highlights that the measure of the traction parameter may be quite experiment-specific, depending on the particular configuration of the gel and surrounding fluid. This is supported by Fig. 2d, where we see the balance between cell traction τ 0 and mixing energy χ that maintains the same equilibrium value of polymer; increasing χ indicates that the gel can equilibrate with a smaller value of τ 0 , and vice versa.
We remark that while we have chosen a particular form for the cell force function G here, other modelling choices have been used. Green et al. (2013), for example, numerically investigated numerous different cell force functions: the Murray force function similar to that described in Sect. 2.2; a 'preferred ECM density' function; and functions incorporating chemical concentration. Green et al. also considered whether contact inhibition should be incorporated in the form λθ p 2 as opposed to the form λn 2 , i.e. acting on the polymer network instead of the cell density. An extension to this work would be to consider different forms of the cell force function, and again, given consistent experimental data to fit the model, to determine if better agreement is found with particular choices.
In the Moon and Tranquillo (1993) study, the gel is constructed as a microsphere; however, in cell-gel experiments, the gels are often thin layers, where the height of the gel is small relative to its length or radius. Therefore, another extension to this model is to study the gel as a two-dimensional thin film, exploiting the ratio of the film's vertical and horizontal length scales as a small parameter to rescale the system and derive a reduced system of equations. In this way, gel behaviour in different experimental settings might be compared and analysed. We aim to present such a model in a future publication.