A two-phase thin-film model for cell-induced gel contraction incorporating osmotic effects

We present a mathematical model of an experiment in which cells are cultured within a gel, which in turn floats freely within a liquid nutrient medium. Traction forces exerted by the cells on the gel cause it to contract over time, giving a measure of the strength of these forces. Building upon our previous work (Reoch et al. in J Math Biol 84(5):31, 2022), we exploit the fact that the gels used frequently have a thin geometry to obtain a reduced model for the behaviour of a thin, two-dimensional cell-seeded gel. We find that steady-state solutions of the reduced model require the cell density and volume fraction of polymer in the gel to be spatially uniform, while the gel height may vary spatially. If we further assume that all three of these variables are initially spatially uniform, this continues for all time and the thin film model can be further reduced to solving a single, non-linear ODE for gel height as a function of time. The thin film model is further investigated for both spatially-uniform and varying initial conditions, using a combination of analytical techniques and numerical simulations. We show that a number of qualitatively different behaviours are possible, depending on the composition of the gel (i.e., the chemical potentials) and the strength of the cell traction forces. However, unlike in the earlier one-dimensional model, we do not observe cases where the gel oscillates between swelling and contraction. For the case of initially uniform cell and gel density, our model predicts that the relative change in the gels’ height and length are equal, which justifies an assumption previously used in the work of Stevenson et al. (Biophys J 99(1):19–28, 2010). Conversely, however, even for non-uniform initial conditions, we do not observe cases where the length of the gel changes whilst its height remains constant, which have been reported in another model of osmotic swelling by Trinschek et al. (AIMS Mater Sci 3(3):1138–1159, 2016; Phys Rev Lett 119:078003, 2017).


Introduction
Tissues in vivo consist of cells living within an extracellular matrix (ECM) comprising a complex network of proteins, glycoproteins, and polysaccharides that surrounds and supports cells within tissues.The ECM provides structural support, forming a scaffold that holds cells together and provides a framework for tissue organisation.It also helps to regulate cell behaviour and tissue function.For example, the ECM contains adhesive proteins, such as fibronectin and laminin, which allow cells to attach to the matrix.This adhesion is essential for cell anchoring and migration, influencing processes like wound healing and tissue regeneration.The ECM acts as a reservoir for growth factors, cytokines, and other signaling molecules.These molecules are stored within the matrix and can be released in response to specific stimuli, influencing nearby cells.It also contributes to the mechanical properties of tissues, such as elasticity, stiffness, and tensile strength (Rozario and DeSimone 2010;Frantz et al. 2010;Humphrey et al. 2014;Dolega et al. 2021).During tissue development and morphogenesis, the ECM guides cell migration, and tissue architecture and shape.It provides spatial cues that direct cell movements and influences the formation of complex tissue structures (Dyson et al. 2016).Importantly, the interactions between cells and the ECM are reciprocal: the ECM is dynamically regulated and can undergo remodeling and degradation by the cells e.g., through the secretion of enzymes called matrix metalloproteinases (MMPs).This process is crucial for tissue repair, regeneration, and remodeling after injury or during normal physiological processes.
In order to recreate a more in vivo-like cellular environment in vitro, cells are often cultured within a gel that emulates (to some degree) the ECM.Collagen, a structural protein that constitutes a major component of the ECM in many animal tissues, is commonly utilized in laboratory studies for this purpose.However, a diverse range of other natural (e.g., Matrigel) or synthetic (e.g., poly(lactic acid)) gels are also employed (Wade and Burdick 2012).Improved understanding of the mechanics of such gels, and of the interactions between cells and the gel, will yield deeper insight into tissue development and functionality.One well-known experiment used to study and quantify how cells remodel the extracellular matrix is the collagen gel contraction assay (Moon and Tranquillo 1993;Barocas et al. 1995).In this assay, cells are cultured within a gel made of collagen, which in turn floats within a bath of nutrient medium.As the cells exert force, they compress the collagen fibres, leading to a reduction in the size of the gel over a period of hours or days.Whilst other gel geometries have been used (e.g., microspheres (Moon and Tranquillo 1993)), a thin disc shape is common, due to ease of fabrication (Vernon and Gooden 2002;Stevenson et al. 2010).Whilst it is expected that cell-induced compaction leads to a reduction in both the diameter and height of the collagen disc, typically the extent of contraction is assessed by measuring only the diameter or the area of the disc (Vernon and Gooden 2002;Stevenson et al. 2010).
In an earlier paper, we developed a model for cell-induced gel compaction incorporating osmotically-driven movement of solvent in or out of the gel (Reoch et al. 2022).This model treats the gel as a two-phase fluid, consisting of polymer and solution, based on previous work (Keener et al. 2011b, a;Mori et al. 2013).However, it additionally includes an equation for the cell density, with the cells exerting a body force on the polymer; mechanical changes in the gel due to this body force were assumed to occur on a short timescale compared with that of cell proliferation and death, so that proliferation and death were neglected.This is reasonable where mechanical changes occur over hours and cell proliferation and death occurs over days.The resulting system of equations was then investigated in a one-dimensional geometry (Reoch et al. 2022).Motivated by the thin geometry commonly used in experiments, in this paper we investigate the behaviour of two-dimensional thin films of gel floating freely in a bath of solvent.We study the 2D Cartesian coordinate system for simplicity of modelling and to facilitate comparison with the 1D Cartesian model presented in earlier work.We continue to neglect cell proliferation and death, assuming the timescale of these process to be much longer than remodelling of the gel under the influence of cell forces.Our two aims are: firstly, to understand the mechanics of the gel in the thin film geometry; and secondly, to compare the emergent behaviours in this system to those in the 1D case.Once more, we hypothesise that the balance between chemical potentials and cell traction stress will be crucial in determining the equilibrium outcome for the gel.To our knowledge, such a model considering the interacting effects of cell traction and osmotic pressure in the thin film geometry has not previously been presented.
Thin film complex fluid models have been employed in a range of biological applications, including modelling cell crawling (Oliver et al. 2005;King and Oliver 2005)), pattern formation in cell aggregation (Green et al. 2017), and osmotically-driven biofilm growth (Trinschek et al. 2016(Trinschek et al. , 2017)).In each case, as in the prototype problem of stretching a thin film of viscous fluid studied by Howell (1996), the small aspect ratio of the film (the ratio of its vertical to horizontal length scales) is exploited to reduce the two-dimensional mass and momentum balance equations to a one-dimensional system at leading order.However, we note that, compared to single-phase fluid models, the two-phase fluid models often require additional assumptions about the relative sizes of some of the parameters in the problem, in order to obtain a one-dimensional reduction.For example, in Green et al. (2017) the parameters describing drag and cell chemotactic effects relative to viscosity were assumed to be large.The chemotaxis scaling encodes the assumption that this is the main driver of cell movement, while the drag scaling couples together the movement of cell and culture medium phases.These scalings were crucial for a non-trivial leading order model to be derived.We make similar assumptions in the model reduction presented here.
This paper is organised as follows.In Sect. 2 we present and nondimensionalise the two-dimensional governing equations and boundary conditions of our model.By exploiting the thin geometry of the film, in Sect. 3 we reduce this model to a system of one-dimensional equations for the leading-order polymer volume fraction and velocity, cell density, film height and film length.Conditions for steady state solutions are then presented in Sect. 4. In Sect.5, we show that when the initial conditions are spatially Fig. 1 Thin film domain = g ∪ s .g is the gel region with θ p > 0, θ s > 0 and cell density n ≥ 0.
s is the region of pure solvent surrounding the gel wherein θ p = n = 0 and θ s = 1.The gel is symmetric about the x-axis and y-axis, and the ratio of gel height to length is small uniform, the thin film model can be further reduced to a single ordinary differential equation for the film height as a function of time, and present numerical solutions in this case.In Sect.7, we consider the behaviour of small spatial perturbations to the initial conditions, before undertaking numerical simulations of the full spatiallyvarying thin film model in Sect.8.The paper concludes in Sect.9 with a discussion of our main results, and suggestions for future work.

Mathematical model
We consider a thin layer of gel consisting of polymer and solvent, which is seeded with cells, and floats freely within a bath of solvent.The set-up is illustrated in Fig. 1.The gel domain is denoted by g , and the surrounding pure solvent domain by s .We adopt a two-dimensional, Cartesian geometry, with spatial coordinates x = (x, y).The gel is assumed to be symmetric about y = 0 and x = 0, with free boundaries at y = ±h(x, t) such that the gel height is 2h(x, t), as well as vertical free boundaries at each end x = ±L(t) giving a gel length of 2L(t).The centre-line of the gel is fixed at y = 0.This allows us to restrict our attention to the region 0 ≤ y ≤ h, 0 ≤ x ≤ L. We assume that all quantities are continuous and differentiable across x = 0 and y = 0. Terms in the external solvent region s will henceforth be denoted with the superscript e; all other terms refer to quantities in the gel region g .The thin gel layer is characterised by a small aspect ratio ε, that is, a small vertical length scale relative to the horizontal length scale (a fact we will exploit in Sect.3).
Our mathematical model is the same as presented in Reoch et al. (2022) and Reoch (2020), though we briefly recapitulate the equations here for completeness.We let the volume fractions of the polymer and solvent be θ p and θ s , respectively, and assume that their mass densities are equal.The cells are assumed to occupy negligible volume, and their density (number per unit volume) is denoted by n. (As explained in Reoch et al. (2022), this latter assumption is based on information from the experimental papers by Iordan et al. (2010) and Moon and Tranquillo (1993) which suggest the volume fraction of the cells when seeded is often O(0.01) or lower.)We assume there are no voids within the gel, so that (2.1) The velocities of the polymer and solvent are denoted by respectively, where v p is the polymer velocity in the horizontal direction and w p is the polymer velocity in the vertical direction, and similarly for the solvent.Cells are assumed to move by advection with the polymer, and by random motion with diffusion coefficient D. Cell proliferation and death are taken to be negligible.Conservation of mass for the polymer, solvent and cells then gives: ) (2.3c) Adding equations (2.3a) and (2.3b) and using the no-voids condition θ p + θ s = 1 gives We model both the polymer and the solvent phases as fluids with a common pressure, P. The polymer is treated as a viscous fluid, where the viscous stress is encapsulated by the stress tensor σ p , related to the rate of strain tensor e p by where the constants η i and κ i , i = p, s, are the dynamic and bulk viscosities, and I is the identity tensor.Whilst similar assumptions are made for the solvent in Reoch et al. (2022), here, for simplicity, we take the solvent viscosities η s and κ s to be zero as in Green et al. (2017), which implies that the solvent stress tensor σ s = 0.As in Keener et al. (2011b) and Reoch et al. (2022), 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), chemical potential gradients and cell-generated forces (which are exerted only on the polymer).We denote the (constant) drag coefficient by ξ , the chemical potentials for the polymer and solvent by μ p (θ p ) and μ s (θ p ), respectively, and the cell force potential by G(n).Conservation of momentum for the polymer and solvent phases then gives (2.6) For the chemical potentials, we use the same forms as Keener et al. (2011b) and Reoch et al. (2022).These are: where f (θ p ) is the free energy per unit volume of gel.In turn, the free energy function is given by 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 and Colby 2003).
Likewise, following Reoch et al. (2022), for the cell potential we set (2.9) The form of this function is similar to that used in earlier works, such as Murray (2001), Moon and Tranquillo (1993), but differs by having a factor of n 2 rather than n in the numerator.This ensures that ∂G/∂n > 0 for all n, which implies that the cell traction force, ∇G = G (n)∇n, is directed up gradients of cell density.The parameter τ 0 ≥ 0 provides a measure for the strength of cell traction forces, and λ ≥ 0 is a contact inhibition parameter, which reduces the force exerted by cells as they become more densely packed.
Taking the sum of equations (2.6) and (2.5) and using the relation θ p ∇μ p = −θ s ∇μ s (which can be derived from (2.7)) we can express equation (2.5) in the form (2.10) Equations (2.6) and (2.10) are used as our momentum balance equations hereafter.

Boundary and initial conditions
Our model consists of equations (2.1), (2.3), (2.6) and (2.10), which require suitable boundary and initial conditions.These will now be specified.Firstly, we consider the motion of the free boundaries at y = h(x, t) and x = L(t), which is given by the standard kinematic conditions (2.12) These conditions imply that polymer 'particles' on the gel's surface must always remain there.We further assume continuity of stress across these free boundaries, hence where the square bracket notation indicates the jump across the boundary.Note that our assumption that the viscosity of the solvent is negligible implies zero tangential stress.Following Mori et al. (2013) we assume that the solvent flux across the gel boundary satisfies where, as before noted, the superscript e indicates a quantity in the solvent domain 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 dependent on the resistance R ≥ 0 of the boundary (an increase in R implies the boundary is less permeable to solvent flow).
Evaluating the force balance equation (2.6) in the solvent domain s , we find that the external pressure P e is at most a function of time, i.e.P e = P e (t); we also note that μ e s = f (0), where f (0) is constant.For the cells, the assumed symmetry of the system about x = 0 implies ∂n ∂ x = 0 at x = 0. (2.15) We impose no-flux cell boundary conditions on y = h and x = L, which implies where n is the unit outward-pointing normal.We assume that the centre of the gel at (0, 0) is stationary.The assumed symmetries of the domain about x = 0 and y = 0 imply that v p must be an odd function of x and an even function of y.They similarly imply that w p is an even function of x and odd function of y.We assume that all the velocities are differentiable throughout the domain.Hence we have: (2.17) (2.18) Note that these conditions similarly hold for v s and w s and their derivatives, although they will not be required for solving the model equations.Initial conditions are required for θ p , n, h, and L. We note that those for θ p , n and h must satisfy the assumed symmetries of the gel about x = 0 and y = 0.A variety of suitable functional forms are considered in subsequent sections.A table summarising the notation used in this paper is given in Appendix A.

Nondimensionalisation and scaling
We now re-scale and nondimensionalise our system to facilitate simplification and analysis.We let L = L(0) be our length scale, N be a characteristic cell density which we choose to be the average cell density at t = 0, and set the time scale to be the ratio of polymer viscosity η p to the free energy scale k B T /ν m .We let H be the height scale, and set this as the average height at t = 0. Thus, the aspect ratio is defined as ε = H/L.We then nondimensionalise the independent variables as follows (where tildes denote dimensionless quantities) with the dependent variables and functions being nondimensionalised thus: The forms of equations (2.1), (2.3a) and (2.3b) are unchanged upon re-scaling, whilst equation (2.3c) becomes where we have introduced the dimensionless diffusion coefficient, (2.20) The momentum balance equations (2.6) and (2.10) (in the x and y directions respectively) now take the form where we have defined the dimensionless parameters (2.23) Note that ξ and τ0 are taken to be O(1).This reflects our assumptions that the unscaled drag and cell traction parameters are large.We now likewise nondimensionalise the boundary conditions.We begin by noting that the kinematic boundary conditions (2.11) and (2.12) are unchanged in form after being re-scaled.On y = h, the stress conditions (2.13) (in the normal and tangential directions, respectively) and the solvent flux condition (2.14) now become: (2.26) where the non-dimensional resistance parameter is defined as and we have introduced the notation P = P − P e and μ s = μ s − μ e s .As with the drag and cell traction parameters ( ξ and τ0 respectively) above, the scaled resistance parameter is taken to be O(1), indicating that resistance is significant.
In order to simplify our notation, we henceforth drop tildes from dimensionless quantities.

Thin film approximation
We now assume that ε 1 and exploit this fact to obtain a simplified, one-dimensional version of our model equations.We begin by expanding our variables in powers of ε: and similarly for θ s , v p , v s , w p , w s , P and L. We then substitute these expansions into our model equations.Unlike previous work discussed in Sect. 1, we expand our variables here in powers of ε instead of ε 2 .This is due to the O(ε) term which appears in interface condition (2.30) after the model is re-scaled.As we will see, the O(ε) terms do not contribute to the leading order model we derive below, but they are included here for completeness.

The y-independence of the leading order dependent variables
We begin by showing that the leading order cell density n 0 , pressure P 0 , polymer fraction θ p 0 and polymer axial velocity v p 0 do not depend on y.This facilitates the derivation of simplified mass and momentum equations, as we will demonstrate subsequently.
Taking equation ( 2.19) at O(ε −2 ), we find At leading order, the no-flux boundary condition (2.31) gives D∂n 0 /∂ y = 0 at y = h 0 .Integrating (3.1) and applying this no-flux boundary condition, we have and so n 0 = n 0 (x, t).(Note that we have assumed D = 0 to derive this result.)We remark that, for our thin-film model to be consistent, this result implies we can only impose initial conditions for n which are independent of y.
where we have defined 0 From interface condition (2.25), at y = h 0 , We have found that P e 0 (t) and higher order equivalent terms do not feature in the final system of equations; therefore, without loss of generality, we now set P e = 0, and as such, P = P for each order of ε.Accordingly, applying (3.5), we find We note that this result holds for all y and, accordingly, throughout the gel.
At O(1), equation (2.21b) yields where μ s 0 = μ s (θ p 0 ); therefore, Now, for (3.8) to hold, we must have that ∂θ Given that G 0 is independent of y and that μ s 0 = μ s (θ p 0 ), both these conditions imply that we must have i.e. at leading order, θ p 0 is independent of y, and accordingly, μ s 0 is as well.Similarly, we have now found that P 0 = P 0 (x, t) by equation (3.6).As noted for n 0 above, for equation (3.9) to be true at all times, it must also hold for the initial polymer fraction θ I (x, y), i.e. our initial condition must satisfy ∂θ I (x, y)/∂ y = 0. Using our definition for 0 together with equation (2.22a) at leading order, we have Noting from (3.6) that 0 = 0, (3.10) becomes and thus F 1 (x, t) = 0. Accordingly, i.e. our leading order polymer axial velocity is independent of y.
From equation (2.21a) at leading order, Therefore, we also have leading order solvent axial velocity v s 0 = v s 0 (x, t).
We now follow the same process to show, in Appendix B, that the O(ε) terms in cell density n 1 , pressure P 1 , polymer fraction θ p 1 , and polymer axial velocity v p 1 are also independent of y.This simplifies later steps in our derivation of a leading-order model.

Derivation of thin film mass balance equations
Having now established that θ p 0 and v p 0 are independent of y, the mass conservation equation (2.3a) at leading order can be integrated with respect to y to give: (where we have used the fact that w p 0 = 0 at y = 0 to determine the arbitrary function of x and t arising from the integration).The kinematic boundary condition at y = h 0 states that Setting y = h 0 in equation (3.15), and using (3.16) then yields This equation describes the mass conservation of polymer over the depth of the thin film.
Similarly, integrating the solvent conservation of mass equation (2.3b) from y = 0 to y = h 0 , and noting that w s 0 | y=0 = 0 gives From interface condition (2.26) at leading order, we also have Using (3.19) and the kinematic boundary condition (3.16) for w p 0 at y = h 0 , we obtain where we have substituted for v s 0 using equation (3.14).Equation (3.20) describes the advection of solvent within the gel.By taking linear combinations of equations (3.20) and (3.17), we can obtain the following equations for h 0 and θ p 0 : In order to derive an expression for n 0 , we evaluate equation ( 2.19) at O( 1).This gives which we express in the following form, noting that all terms on the right-hand side are independent of y: After integrating between y = 0 and y = h 0 , we have We note that ∂n 2 /∂ y = 0 at y = 0 from the symmetry condition (2.15).From the no-flux boundary condition (2.31) at y = h 0 , we find that at O(ε 2 ), On substituting in these expressions and rearranging, equation (3.25) becomes Given that h 0 > 0 in the gel, we divide by h 0 and use (3.15) to substitute for w p 0 , obtaining the mass conservation equation for n 0 ,

Derivation of an equation for v p 0
Having derived equations for the polymer volume fraction, cell density and the height of the film, we now close our system of leading-order thin-film equations by finding an expression for Next, we integrate equation (3.29) with respect to y from y = 0 to y = h 0 to obtain (3.31) Further, by symmetry at y = 0, we have w p 0 (x, 0, t) = 0, and so ∂w p 0 /∂ x = 0 at y = 0. We also have that v p 2 is an even function of y about y = 0, and thus ∂v p 2 /∂ y = 0 at y = 0. Therefore, from equation (3.30), Applying Leibniz' integral rule in equation (3.32), after a little algebra we obtain 123 Our next step is to eliminate the 2 term in (3.33).To do this, we consider (2.22b) at O(1), finding where F 3 (x, t) is an arbitrary function.At O(1), boundary condition (2.25) on y = h 0 states that Therefore, we find F 3 (x, t) = 0, and from (3.34), We can now simplify the integral term in (3.33).After substituting in the expression for 2 above and noting that the resulting integrand is independent of y, equation (3.33) becomes Now, using equation (3.15) to substitute for the w p 0 term, and integrating with respect to x, we find Since the left-hand side of this expression is independent of y, we must have F 4 = F 4 (t).On substituting from equations (3.15) and (3.36) into equation (3.38) we obtain The interface condition (2.28) at x = L 0 then supplies Thus, we find F 4 = 0, and equation (3.38) simplifies to which gives a leading order expression for the polymer axial velocity v p 0 .We have thus obtained a closed system of equations for h 0 , θ p 0 , n 0 and v p 0 -namely, (3.21), (3.22), (3.28) and (3.41) -which constitutes our simplified thin-film model.Finally, we note from the leading order interface condition (2.30) at x = L 0 that v s 0 = v p 0 .Substituting this into (3.14),we find, at x = L 0 , We also have no cell flux at x = L 0 , therefore ∂n 0 /∂ x = 0, and we can express equation (3.42) as this indicates that at x = L 0 we must have ∂θ p 0 /∂ x = 0.For consistency, we will always use initial conditions such that this boundary condition is satisfied.We note that physically this may not always be true; this may lead to more complicated scenarios (such as the existence of boundary layers) which we do not address here.

Summary of thin film model equations
Given the length of the calculations in the preceding section, for convenience we now summarise the system of thin film model equations which we consider in the remaining sections of this paper.This comprises (3.22) for θ p , (3.21) for h, (3.28) for n and (3.41) for v p .(Note that henceforth, we drop the zero subscript denoting leading order quantities to reduce notational clutter.)We eliminate ∂θ p 0 /∂t from equation (3.41) using (3.22), and after a little algebra, our system becomes: The above system is subject to the boundary conditions and initial conditions Note that we take the initial condition θ I (x) to be differentiable and to satisfy ∂θ I (0)/∂ x = 0, so that for all time, ∂θ p /∂ x| x=0 = 0. We similarly take the initial height h I (x) to be differentiable and satisfy ∂h I /∂ x| x=0 = 0, so that ∂h/∂ x| x=0 = 0 for all time (required by the symmetry of the problem).Finally, we note that equation (3.43) implies ∂θ p /∂ x| x=L(t) = 0, so our initial conditions must satisfy:

Steady state conditions
We now consider the steady state solutions of our thin film model (3.44).We being by returning briefly to the conservative form of the polymer advection equation (3.17).At equilibrium, this supplies where we have integrated with respect to x and applied the zero velocity boundary condition at x = 0. Given θ p > 0 and h > 0 within the gel, this means that we must have v p = 0 at equilibrium.Since L = v p | x=L(t) , we have L = L * (a constant).Then, from equation (3.14), we have We note that since v p = 0, and v s = 0 at x = 0 by symmetry, we require v s = 0 throughout the gel, as there will otherwise be motion of the solution relative to the polymer, and θ p will vary in time.Hence, the right-hand side of (4.2) must be zero at a steady state.Using these pieces of information, equation (6.1d) implies that we must have We now turn to the cell advection-diffusion equation (6.1c), which reduces to Assuming D > 0, integrating (4.4) and applying the no-flux boundary condition (6.2b), we find, given h > 0, this indicates that, at equilibrium, and accordingly, n must be spatially uniform.
Using the fact that G = G(n), equation ( 4.3) now gives and therefore, we must have Therefore, necessary and sufficient conditions for equilibrium in the thin film are spatially uniform values θ p = θ * and n = n * that satisfy (4.3).We note that there are no restrictions on h, i.e. the height can be non-uniform in space at the gel's steady state.The equilibrium conditions (4.3) match those for a one-dimensional gel with D = 0, as presented in (Reoch et al. 2022).Therefore, any equilibrium solution for θ p and n in the 1D Cartesian model with D > 0, is an equilibrium of the thin film model.

Reduced model for uniform initial conditions
Given spatially uniform initial conditions for θ p , n and h, we now show that we can simplify the thin film system of equations (3.44) and boundary conditions (3.45) to an even simpler form.In this case, our model reduces to an ODE for h as a function of time, with the other variables specified in terms of h.We begin by setting the initial conditions to be where we have scaled h and n on the initial height and cell density respectively.Noting that for uniform initial conditions, the free energy and cell force are also initially uniform in x, this means that the velocity equation (6.1d) at initial time is Given that the terms on the right hand side are all initially independent of x, the initial velocity is where we have used the condition v p = 0 at x = 0. We can then verify from equations (3.44a), (3.44b) and (3.44c) that θ p , h and n are initially independent of x.By making the ansatz that v p is linear in x, and θ p , h and n are functions of time only, on substituting into (3.44a),(3.44b) and (3.44c) we arrive at the system of ODEs: ) ) (5.4d) The kinematic condition gives (5.5) On applying the initial conditions, we can find θ p , v p , n and L in terms of h as: (5.6d) Having found solutions θ p and n as functions of h, we can also express the chemical potential μ s (θ p ) and cell force function G(n) as functions of h, and accordingly, reduce our model to the following single ODE for h(T ), From this reduced model, we can draw a number of conclusions about how the gel behaves under spatially uniform initial conditions.Firstly, that if there is no initial spatial variation in the volume fractions, cell density and gel height, the gel will remain uniform in space for all time.We see that the length of the gel is equal to its height as it evolves, therefore any time-variation in the free boundaries of the gel will occur in the same manner both length-wise and height-wise.We also see that the polymer fraction and cell density are inversely proportional to the squared height, scaled by the initial conditions.Finally, the reduced model is independent of drag; with no dependence on the x-coordinate, there is no relative motion between the polymer and solvent as the gel evolves and so no shearing takes place.
The entire system is driven by the balance between the solvent chemical potential μ s and the cell potential θ p G inside the gel with the external chemical potential μ e s .If these forces are in balance, the system is in equilibrium, as expected.For the height to be increasing, and accordingly, the gel to be swelling, we require that μ e s is greater than the sum of θ p G and μ s ; the opposite holds for the gel to contract.The gel equilibrates when these forces are in balance, that is the same condition as derived in Sect. 4 above.The rate at which the gel evolves is determined by equation (5.7).Accordingly, increasing the resistance of the interface R slows the rate of change of h and the rest of the system.Similarly, larger values of parameters such as the mixing energy χ and cell traction τ 0 appearing in μ s and G respectively will increase the rate of gel evolution.

Numerical simulations of the reduced model
We solve the reduced model described by equations (5.6) and (5.7) using the inbuilt MATLAB ODE solver ode15s which is designed to handle stiff systems.For a gel without cells, we can see either swelling or contraction take place, depending on the chemical potentials in the system.Figure 2a demonstrates swelling induced by osmotic pressure (with mixing parameter χ = 0.75), while, conversely, Fig. 2b shows a case where the gel contracts (with mixing increased to χ = 1.5, promoting separation between the polymer and solvent).The equilibria reached here, θ * = 0.45 for the swelling case and θ * = 0.86 for the contracting case, are, of course, the same as found for the equivalent initial conditions and parameters in the 1D Cartesian version of the model (Reoch et al. 2022;Reoch 2020).
Introducing cells into this system can precipitate a switch to contraction in a gel that would otherwise swell.Note that this is due to the traction forces exerted by the cells; the spatial uniformity of the cell population means that diffusion plays no role in the evolution of the system.We take the same parameter values as in the swelling gel above (Fig. 2a) and introduce a cell population (n I = 1, τ 0 = 1).In Fig. 3a, we see that this gel now contracts due to the presence of cells, reaching a steady state with a significant reduction in height and length.This highlights that in this thin film geometry, given sufficient traction stress (τ 0 = 1), the presence of cells can outweigh the osmotic swelling pressure created by chemical potentials in a gel.Reducing the traction parameter to τ 0 = 0.1 (Fig. 3b), we see that the presence of cells does not necessitate contraction; rather it is evident that, due to a weak cell contribution, osmotic pressure is still the dominant driver of the gel's behaviour, and the thin film still swells Fig. 2 Time evolution of a cell-free thin gel with uniform initial conditions; θ p is shown in blue, and h in maroon, with h = L. Common parameters θ I = 0.6, n I = 0, h I = 1, N = 100, R = 1.With χ = 0.75 (Fig. 2a) the gel swells uniformly across its domain to the equilibrium (θ * , h * , L * ) = (0.45, 1.16, 1.16).With χ = 1.5 (Fig. 2b) the gel contracts to the equilibrium (θ * , h * , L * ) = (0.86, 0.84, 0.84) to an equilibrium.The equilibrium values of polymer and cell density again match those found for the 1D case for the same initial conditions and parameters (Reoch et al. 2022;Reoch 2020).
Increasing the resistance parameter R slows the evolution of the gel and hence increases the time taken to reach a steady state; however, it does not affect the eventual equilibrium reached.Figure 3c shows the effect of increasing the resistance parameter from R = 1 to R = 5, which increases the time taken for the gel to reach its steady state approximately five-fold, from T ≈ 2 to T ≈ 10.This relation between R and T is expected, given that the solution for h in equation (5.7) is of the form F(T /R).
For both the contracting and expanding examples here, the gel remains spatially uniform throughout its evolution to steady state.The height and length are necessarily equal by (5.6c), indicating that the gel grows in a uniform ratio horizontally and vertically.Next, we will explore whether this remains the case when the gel is constructed with non-uniform initial conditions.

Transformation to a 1D fixed domain
In order to facilitate the analysis and numerical solution of our model in spatiallyvarying cases, we transform to a fixed domain using the coordinate transformation t = T , x = L 0 (T )X .The equations (3.44) then become: Fig. 3 Time evolution of a thin cell-gel system with uniform initial conditions.In Fig. 3a and b, θ p is shown by the solid blue curve, h = L is the solid maroon line, and n is the dotted purple line, for common parameter values θ I = 0.6, n I = 1, h I = 1, χ = 0.75, N = 100, R = 1, λ = 1.With τ 0 = 1 (Fig. 3a) the gel contracts to the equilibrium (θ * , n * , h * , L * ) = (0.86, 1.44, 0.83, 0.83) due to the presence of cells.With τ 0 = 0.1 (Fig. 3b) the gel swells to a steady state due to osmoticpressure counteracting weak cell traction, (θ * , n * , h * , L * ) = (0.54, 0.91, 1.05, 1.05).Figure 3c compares the effect of the size of the resistance parameter on the evolution of the gel height; R = 1 is shown by the light blue dotted line, R = 5 by the maroon dashed line, parameter values otherwise as for Fig. 3a 2h L subject to boundary conditions and initial conditions Note that our assumptions regarding the symmetries, etc. of the initial conditions (3.47) now become:

Small time evolution of spatially perturbed uniform equilibria
We begin our investigation of the behaviour of the model with non-uniform initial conditions by considering the short time behaviour when equilibrium initial conditions are subject to spatial perturbations.This analysis will provide insight into the stability of equilibria: unstable equilibria will see an increase in the amplitude of the perturbations, while for stable equilibria, the amplitude of the perturbations will decay.As seen in Sect.4, all equilibria have spatially uniform polymer fraction (θ p ) and cell density (n) but they need not have spatially uniform film height (h).However, in this section we restrict our attention to equilibria with spatially uniform h, to facilitate finding analytic solutions and, hence, simplify our analysis.We denote the dimensionless equilibrium values by asterisks, L * , θ * , n * , h * , v * (where v * = 0).The equilibrium values of h, n and L are used as the characteristic values to scale these variables; this means that h * = n * = L * = 1.We introduce a short timescale, T , such that T = δ T where δ 1, and expand our solutions as power series in δ: with expansions for θ p , h and n similar to that for v p .We take the spatial perturbations to have an amplitude , where δ 1. (Note that the amplitude is distinct from the aspect ratio ε exploited in Sect.3.) We then take the series (7.1), etc., and expand each of the terms L j , v j , θ j , h j n j , j = 1, 2, . . . in powers of , for example, We take the initial conditions to be Note that we set L 0 = L * = 1, i.e. we do not perturb the initial length of the gel from its equilibrium value.Since our initial conditions are equilibria, we have v 00 = 0, with higher order terms of v 0 determined through analysis of the momentum balance equation (6.1d).We set θ 01 = cos(γ X ), n 01 = N 01 cos(γ X ), and take h 01 = −H 01 cos(γ X ), where N 01 and H 01 are constants that are O( 1).(Given that we expect h to change in the opposite way to θ p and n, we have taken h 01 to have the opposite sign to that of n 01 and θ 01 without loss of generality.)As required, θ 0 , n 0 and h 0 satisfy the symmetry boundary conditions (6.2a), (6.2b) and (6.2c) at X = 0 for any choice of γ .The no-flux boundary conditions (6.2a) and (6.2b) at X = 1 require that γ = Z π for some positive integer Z .Finally, we wish to ensure that, for any choices of the constants H 01 and N 01 , the masses of polymer and cells under the perturbed initial conditions are, to O( ), equal to the masses for the unperturbed initial conditions (which are θ * for the polymer and 1 for the cell density).Integrating θ 0 h 0 L 0 across the spatial domain 0 ≤ X ≤ 1, we have Since sin(γ ) = 0 for all valid choices of γ , we see that mass is conserved to O( ), regardless of our choice of H 01 .On evaluating the integral of n 0 h 0 L 0 , we find that we are similarly free to set N 01 to any O(1) value.
To find v 01 , we use equation (6.1d) at O( ), obtaining the expression After substituting for θ 01 and n 01 , this simplifies to where We therefore have the solution Solving the mass conservation equations (6.1b), (6.1c) at O(δ), we find that θ 10 = n 10 = h 10 = 0, as these terms depend on v 00 = 0.The kinematic boundary condition (6.2e) for L similarly gives L 10 = 0. We therefore consider the mass conservation equations at O(δ ).From equation (6.1b), we find Similarly, we find from equations (6.1a) and (6.1c) respectively that Finally, from the kinematic boundary condition at O(δ ), we find that L 11 = v 01 | X =1 = 0. Thus, we have the small time analytic solutions (7.7d) We note that these solutions satisfy the no-flux boundary conditions at X = 1.The behaviour of the amplitudes of the perturbations can be determined by considering the square-bracketed terms in equations (7.7).We see that the growth or decay of the perturbations of θ p and h is governed by the sign of z.For z > 0, the magnitudes of the spatial perturbations are decreasing in time, and accordingly, the local spatial variations decay.For z < 0 on the other hand, the amplitudes of the polymer fraction and height perturbations grow, indicating that the equilibrium is unstable.The cell density (7.7b) follows this same behaviour, although the presence of diffusion can be stabilising, depending on the balance of parameters.The stability conditions for spatially uniform equilibria of the thin film model, including the value of z, are thus equivalent to those presented in Reoch et al. (2022) for the one-dimensional case; we refer readers to the stability diagrams there given and do not give similar figures here.

Thin film numerics
We now present our numerical scheme to solve the thin film equations (6.1) in the spatially non-uniform case.We being by re-writing the equations for θ p and h in conservative form by defining the quantities Q(X , T ) = θ p h and W (X , T ) = θ s h = (1 − θ p )h, and noting the following relations: We then transform the chemical potential μ s into a function of Q and W . Equations (6.1a), (6.1b) and (6.1d) can then be expressed respectively as: where The cell advection-diffusion equation remains as given in equation (6.1c).Hence, our numerical code solves the system given by (8.1a) -(8.1c) together with (6.1c).
For this, we implement a finite difference scheme in MATLAB which discretises the system of equations using a uniform spatial grid between X = 0 and X = 1.The force balance equation (8.1c) is first-order in velocity v p , and therefore we use a cumulative trapezoidal scheme to numerically integrate the expression across the spatial domain and update the velocity at each new time step.Central differencing is used for spatial derivatives in (8.1c), except for the derivatives of (8.2a) at the endpoints of the domain, where one-sided differences are used.A Crank-Nicolson method is used to solve equations (8.1a), (8.1b) and (6.1c).The end time is chosen to be large enough that the gel reaches a steady state or θ p approaches 0 or 1 (in which case our model breaks down).
We find that this scheme conserves mass effectively.Using a time step dT between 10 −6 and 10 −5 and spatial step d X = 0.025 in the simulations which follow, the worst-case change in mass between initial time and end time for the cell density or polymer fraction was 2.69 × 10 −5 %.To check our numerical scheme, we compared the full numerical solution with the small time solution detailed in Sect.7, and also compared the solution obtained with this code for uniform initial conditions with the ODE solution from Sect.5.1, finding good agreement in each case (results not shown).While our numerical scheme has performed well over a wide range of parameter choices and initial conditions, we note that some instability and non-convergence has been encountered for particular parameter combinations where the gel evolution is rapid, e.g. with large values of τ 0 or χ .We do not consider such cases here.Throughout the simulations presented in this section, we keep certain parameters fixed and study the effects of changing others between simulations.The fixed parameters and their values are given in Table 1.The values used for the parameters and initial conditions which are varied between examples are presented in Table 2.We note that the initial conditions for θ p , n and h will have spatially varying components included on occasion in the simulations which follow; accordingly, we scale the initial conditions for height and cell density using the average initial value for each, such that the mean values h I = 1 0 h I (X )d X = 1 and nI = 1 0 n I (X )d X = 1 or 0 depending on the presence or absence of cells.Similarly, when spatial perturbations are added to θ p , these are such that θI = 1 0 θ I (X )d X.

Non-uniform initial conditions
We now consider spatially varying initial conditions in one or more of the polymer fraction, cell density and height.We analyse the gel behaviours emerging in different cases, in particular evaluating whether spatial non-uniformities persist or are smoothed out as the gel evolves over time, and whether perturbing different dependent variables results in different impacts on the final outcome.

Spatially varying initial polymer, cell-free gel
We first consider a cell-free gel with a non-uniform initial polymer distribution.We take the initial condition θ = 0.6+0.02cos(π X ), noting that θ I satisfies the boundary conditions ∂θ p /∂ X = 0 at X = 0, 1, and that the initial mass of polymer in the gel remains 0.6 under this initial condition.This resembles a gel where the polymer fraction is initially slightly larger around the centre of the thin film at X = 0. We will firstly discuss the equilibrium outcomes for the gel, then describe how the variables change over time.
We see in this environment that spatially varying steady states can be found without the presence of cells.Figure 4a displays the time evolution of θ p (X , T ) at X = 0, 1, while Fig. 4b shows the time evolution of h(X , T ) at X = 0, 1 as well as L(T ).We see that the gel height h evolves to a non-uniform equilibrium here.The mixing parameter χ = 0.75 is at a value that promotes mixing between the polymer and solvent.This drives an osmotic pressure gradient, causing the gel to swell to an equilibrium state, with solvent flowing into the gel.The polymer fraction converges to a spatially uniform value as it evolves, reaching a steady state where θ * = 0.45.Meanwhile, we see that non-uniformities develop in the gel height; these spatial variations persist at equilibrium where the mean equilibrium value h * = 1.16 and the amplitude A h * = 0.019, where A h * = (h * max − h * min )/2.The amplitude A h * is of a similar magnitude to the amplitude of the initial polymer fraction.The gel length at equilibrium L * = 1.16 is equal to the mean equilibrium value of the height h * .
Figures 4c and d show the spatial distributions of θ p and h respectively across the gel length at increasing points in time.As seen in Fig. 4c, the initial non-uniformity in the polymer distribution quickly smooths out so that θ p is uniform across the spatial domain.Meanwhile, Fig. 4d demonstrates that sinusoidal variations matching the shape of those in the initial polymer arise in the gel height.These variations persist over time, resulting in varying height at the gel's steady state.
In response to the osmotic pressure gradient here driven by the free energy, we see more solvent enter the gel over early time (e.g.T = 0 to T = 0.5) in the regions of higher polymer fractions near X = 0, resulting in these areas of the gel becoming locally thicker.This is seen in the gel height increasing to a greater degree close to X = 0 since more solvent is entering the gel in that region.Conversely, we see θ p increase near X = 1 over this time period, i.e. there is some localised contraction in the gel due to the initial presence of more solvent in this region.This corresponds to decreases in h seen at corresponding times.By T = 1, the gel swells across the spatial domain, with a uniform polymer profile developing as solvent continues to enter the gel more rapidly in areas of greater polymer concentration.The height continues to increase as the gel swells, maintaining its non-uniform distribution.These variations that develop and persist in h correspond to local variations in mass across the spatial domain that exist from the initial non-uniformity in θ I .Accordingly, while the fraction of polymer is constant by the time the gel equilibrates, the mass of polymer per unit length θ * h * (X ) varies in space.
In Sect.4, we found that ∂θ p /∂ X = 0 is a necessary condition for equilibrium in the thin film.We see here that the polymer is redistributed evenly such that there is a balance in chemical potential from gel to solvent and within the gel itself.The extra mass at different points in space (coming from the spatially varying height) enables a constant Fig. 4 Time evolution of a cell-free gel with non-uniform initial polymer fraction θ I = 0.6+0.02cos(π X ), and parameters n I = 0, h I = 1, χ = 0.75, ξ = 1, R = 1, expanding to equilibrium (θ * , h * , L * ) = (0.45, 1.16, 1.16).Figure 4a shows θ p versus time T at X = 0 (solid blue curve) and X = 1 (dashed red curve); the curves are identical to graphical accuracy.Figure 4b shows h versus T at X = 0 (dashed light blue curve) and X = 1 (dotted maroon curve), along with L(T ) as the solid gold curve.Figure 4c and d show curves θ p versus X and h versus X , respectively, at times T = 0, 0.1, 0.2, 0.5, 1, 2, 8, 120, with time increasing in the direction of the black arrow.The polymer fraction evens out to a uniform equilibrium as the gel swells (Fig. 4a and c).Spatial variations develop in the height in response to the initial non-uniform polymer distribution and persist to equilibrium (Fig. 4b and d) polymer fraction to be maintained at equilibrium.We note that this model does not consider surface tension.With surface tension present, we might expect the variations in height to smooth over time as well; however, in its absence, there is no force driving the surface to flatten out and we see the non-uniformities persist at equilibrium.
We see the same qualitative outcome for the gel when taking non-uniform h I with uniform θ I , i.e. at the resulting steady state, h * will vary in space while θ * is uniform (results not shown).As seen in the previous example, the variations in mass across the spatial domain allow for a uniform polymer fraction to be maintained at equilibrium.

Spatially varying initial polymer, cell-gel system
We now take the system presented in Sect.8.1.1 and introduce a cell population where n I = 1 and τ 0 = 1.We maintain the non-uniform initial condition for polymer, θ I = 0.6 + 0.02 cos(π X ). Figure 5 displays this system's evolution.While the gel Fig. 5 Time evolution of a cell-gel system with non-uniform initial polymer fraction θ I = 0.6 + 0.02 cos(π X ), and parameters n I = 1, h I = 1, χ = 0.75, ξ = 1, R = 1, τ 0 = 1, D = 1, contracting to equilibrium (θ * , n * , h * , L * ) = (0.86, 1.44, 0.83, 0.83).Figure 5a shows that spatial variations in the polymer profile decay quickly over time, whilst in Fig. 5b small spatial variations briefly emerge in the cell density, but dissipate before the gel equilibrates.By contrast, in Fig. 5c spatial variations which emerge in the gel height increase in magnitude throughout to equilibrium.Profiles are plotted at T = 0, 0.1, 0.2, 0.5, 0.8, 1.2, 1.6, 3, with time increasing in the direction of the black arrow was previously seen to swell, the introduction of cells switches the gel's behaviour to contraction, with the forces the cells generate outweighing the chemical potential gradient.The gel reaches a steady state where θ * = 0.86, n * = 1.44, h * = 0.83, L * = 0.83.We note that this is the same equilibrium in θ * and n * as obtained using the same parameter values in our one-dimensional model (Reoch et al. 2022).
As in Sect.8.1.1,the non-uniformity in θ p evens out over time, with h developing spatial variations that remain present at equilibrium; Fig. 5a and c show how the spatial profiles of θ p and h respectively change over time.Meanwhile, variations also appear in the cell density n while the gel contracts; the cell density increases more around X = 1 in response to the smaller initial polymer fraction there.However, as time progresses, the variation in n decays due to the presence of diffusion (see Fig. 5b).

Spatially varying initial cell density, cell-gel system
We now take the initial cell density to be spatially varying, such that n I = 1 + 0.02 cos(π X ), with a uniform initial polymer fraction θ I = 0.6.Small spatial variations arise in both the polymer fraction and height here as the gel evolves (see Fig. 6)..86, 1.44, 0.83, 0.83).Spatial variations briefly emerge in the polymer fraction (Fig. 6a) and the gel height (Fig. 6c), which dissipate before the gel equilibrates, whilst spatial variations in the cell density decay quickly over time (Fig. 6b).Profiles are plotted at T = 0, 0.06, 0.1, 0.2, 0.5, 0.8, 1.2, 3, with time increasing in the direction of the black arrow The non-uniform cell distribution, shown in Fig. 6b, leads to greater forces initially being applied in the negative X -direction; this cell traction induces spatial gradients in the polymer profile and, accordingly, the height, as more solvent is forced from the gel.Figure 6a shows θ p increasing towards X = 0 over early time in response to the cell force gradient, while in Fig. 6c, we see a corresponding decrease in height around X = 0. Contrary to the previous examples seen in this Section, the non-uniformity that arises in the height as the gel evolves does not continue to equilibrium.In this instance, with the polymer initially uniform, the small local variations in the gel's thickness do not persist at equilibrium.As required, the cell density is constant when the gel equilibrates, with the strong diffusion coefficient playing a significant role in smoothing out the initial variations.The gel reaches the same steady state as the previous example with θ * p = 0.86, n * = 1.44, h * = L * = 0.83; therefore, we see that varying the initial cell distribution does not led to greater contraction in the gel.

Spatially varying initial height, cell-gel system
We now vary the initial height for this gel, such that h I = 1 + 0.02 cos(π X ), while taking θ I = 0.6 and n I = 1 to be constant.The gel again contracts to an equilibrium with θ p and n being spatially uniform, while spatial variations in h persist through to the gel's steady state.The mean equilibrium values of the model variables are unchanged from previous examples.As this example does not demonstrate any new qualitative outcomes, we do not include any figures.
We note that taking different combinations of spatially varying initial conditions in h, n and θ p will result in the same qualitative outcomes for the system -if θ p or h are spatially varying initially, then h will be spatially dependent at equilibrium, regardless of the initial cell density; only varying n initially will not induce a nonuniform equilibrium by itself.

Influence of drag and resistance
We now investigate the effect that the drag parameter ξ and the resistance parameter R have on the gel's evolution.We take a contracting gel with cells, with parameter values and initial conditions as given in Fig. 5, where the initial polymer fraction is non-uniform.We reiterate that in this case, the gel will equilibrate to a uniform value of polymer, θ * = 0.86.We now modify the drag parameter ξ .In Fig. 7a and b we compare the polymer fraction's spatial evolution over early time (up to T = 0.8) for high drag (ξ = 4) and low drag (ξ = 0.2) respectively.In the high drag case, we see little change over this time to the spatial structure of the gel as it contracts.The amplitude of the variations in the polymer fraction decreases over time (amplitude A θ p = 0.007 at T = 0.8); however, it retains the sinusoidal shape of the initial condition.In this case, the large drag coefficient slows down solvent flow through the gel in the x-direction.It is therefore easier for solvent to flow out of the gel primarily in the thin direction.It does this at a relatively uniform rate across the domain, hence the spatial distribution only slowly decreases in amplitude.With low drag on the other hand, we see that the polymer fraction changes more rapidly near X = 1 than in the high drag case.The polymer moves quickly towards a uniform spatial distribution, since it is now easier for fluid to flow longitudinally within the gel as well as in the vertical direction.With more cell forces being applied in the area near X = 1, the polymer fraction increases at one point (at around T = 0.4) as it evolves, before evening out again as the gel moves towards its equilibrium state.We note that the polymer fraction reaches a uniform steady state at T ≈ 3 in both the low and high drag cases here (results not shown).
The resistance parameter R affects the speed at which fluid can flow across the gel-solvent interface at both y = h and X = 1.We see dramatic differences when comparing the evolution in θ p for low resistance with R = 0.4 in Fig. 8b with high resistance, R = 4 in Fig. 8a (we note that in both these examples ξ = 1).For low resistance, the gel quickly contracts, with θ p smoothing out as it moves towards its steady state.With this small value of R, the gel has equilibrated by T = 0.8.For high resistance, the evolution is significantly slower.Indeed, by T = 0.8, the polymer Fig. 7 The effect of varying the drag parameter, ξ .For large drag (ξ = 4, Fig. 7a) spatial variations in polymer density slowly recede as the gel contracts over time.By contrast, in the case of low drag (ξ = 0.2, Fig. 7b), the spatial variations quickly smooth out as the gel contracts.Profiles are plotted at T = 0, 0.02, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, with time increasing in the direction of the black arrow.Initial conditions and parameter values otherwise as given in Fig. 5 Fig. 8 The effect of varying the resistance parameter, R. When resistance is high (R = 4, Fig. 8a) the polymer fraction decreases very slowly due to the impermeability of the boundary, whilst when the boundary is more permeable (R = 0.4, Fig. 8b), contraction is much faster.Profiles are plotted at T = 0, 0.02, 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, with time increasing in the direction of the black arrow.Initial conditions and parameter values otherwise as given in Fig. 5 fraction has only increased to approximately θ p = 0.64 at its maximum.Interestingly, in this case, while the fraction of polymer is only changing slowly, the amplitude has halved from its initial value by T = 0.8.This indicates that, while the resistance is slowing the flow of solvent across the gel's boundary, solvent inside the gel is flowing quickly enough to flatten the polymer's distribution.

Zero-diffusion case
The derivation of the thin film model, wherein the diffusive flux terms are used to derive that n is independent of y, requires D to be O(1).Accordingly, the examples presented in this chapter so far have been generated with D = 1.Nevertheless, for the sake of completeness, we now investigate the outcomes for the system without diffusion, although we cannot guarantee the validity of our thin film reduction in this case.With zero diffusion, cells must move with the polymer and the cell distribution does not need to be uniform at steady state, indicating that θ p in turn is also not required to be uniform at equilibrium.
By way of example, we take D = 0 with initial conditions n I = 1 + 0.02 cos(π X ), θ I = 0.6, h I = 1.We note here that due to issues with code convergence, we have run these examples with a smaller interaction energy than previously, taking χ = 0.4.We show in Fig. 9a, b that with no diffusion and a non-uniform initial cell density, we obtain equilibria that are non-uniform in the cell density and polymer fraction.This reflects similar examples in our 1D model (Reoch et al. 2022).In this instance, the gel reaches a mean equilibrium cell density n * = 1.29 with amplitude A n * = 0.033, which is greater than the initial amplitude.The polymer fraction evolves to mean equilibrium value θ * = 0.78 with amplitude A θ * = 0.004, while for the height, h * = 0.88 with amplitude A h * = 0.002.Figure 9c demonstrates the time evolution of velocity v p at different points in the spatial domain.We see that the velocity goes to zero at all shown spatial points, confirming that the gel is at a steady state.With only cells taken to be initially non-uniform, the gel reaches a steady state where the cells, polymer and height are all non-uniform.This is a significant difference to the case with diffusion, where non-uniform initial cell profiles did not lead to spatially varying steady states.While this simulation does not fit with the particular derivation of the thin film system here, it does demonstrate the possibility that spatially dependent solutions may occur in both the polymer and cells.Examples with a small diffusion coefficient (e.g.D ≈ 1x10 −3 − 1x10 −5 ) were found to reach a similar non-uniform quasi-steady state; however, given the presence of the small diffusive flux, cells and polymer moved very slowly towards a uniform distribution (results not shown).

Discussion
In this paper, we have presented a model to study the behaviour of a thin film of gel floating freely within a bath of solvent.Starting from the equations of our earlier model of cell-induced gel contraction (Reoch et al. 2022), we exploited the small aspect ratio ε of the problem, assumed symmetry about x = 0 and y = 0, and showed that under a particular set of scalings, the original two-dimensional system of equations can be reduced to a one-dimensional system of four coupled PDEs.This leading-order model consists of equations for the gel half-height h (3.21), polymer fraction θ p (3.22), cell density n (3.28) and polymer velocity v p (3.41).Such a model, considering cell traction stresses and osmotic pressure in this thin film setting has not, to our knowledge, been presented previously.
The thin film model presented here has some key differences to the 1D model presented in Reoch et al. (2022).In the previous geometry, all solvent flowed in or out of the gel horizontally through the gel's endpoints at X = ±1.In contrast, in the thin film model, most of the solvent flow into or out of the gel occurs across the long boundaries at y = ±h.This manifests itself in the additional terms seen in equation (3.20), which describes (at leading order) mass conservation of the solvent Fig. 9 Time evolution of a cell-gel system with zero cell diffusion, a non-uniform initial cell density n I = 1 + 0.02 cos(π X ), and parameters θ I = 0.6, h I = 1, χ = 0.4, ξ = 1, R = 1, τ 0 = 1, D = 0, contracting to equilibrium ( θ * , n * , h * , L * ) = (0.78, 1.29, 0.88, 0.88).Spatial variations in the cell profile grow and persist at equilibrium (Fig. 9a), whilst similar spatial variations in the polymer density also emerge and persist (Fig. 9b); profiles are plotted at T = 0, 0.05, 0.1, 0.2, 0.4, 0.8, 1.6, 8, with time increasing in the direction of the black arrow.Figure 9c shows curves of the velocity v p , versus time T , at X = 0 (blue), X = 0.25 (red), X = 0.5 (yellow), X = 0.75 (purple), X = 1 (green); v p goes to zero across the spatial domain as the gel reaches its spatially varying equilibrium fraction across the thickness of the gel, and includes a source term which depends on the balance of chemical and cell potentials, scaled by the resistance of the interface to fluid flow.This source term effectively emerges as a result of integrating over the height of the gel.Another key difference is that we have an equation (3.41) for the polymer velocity in the thin film model which is first order in space, as opposed to a second order equation previously.Finally, we note that the thin film approximation relies on the presence of the cell diffusion terms to show that n 0 is independent of y; therefore, unlike the previous 1D case, we cannot set D = 0 without violating one of the assumptions underpinning the thin film model.Despite these differences in the detail of the equations, the conditions for a steady state remain the same as in the one-dimensional model of Reoch et al. (2022).
For spatially uniform initial conditions, we have shown that the solution of thin-film equations is also spatially uniform in the dependent variables θ p , n and h.In fact, the model can be reduced to a single ODE in time for the height h.By contrast, in the 1D Cartesian model such spatially-uniform solutions only occur for no drag, ξ = 0 (Keener et al. 2011b;Reoch et al. 2022).This is because, in the 1D case, fluid must flow into or out of the gel through the boundary at X = 1, producing non-uniform spatial profiles in the dependent variables as it moves through the gel (these non-uniformities generally smooth out by the time the gel equilibrates).In the thin film, solvent can flow across the long boundary at y = h and so drag has less influence on the gel's behaviour.The dynamics of the ODE model are driven by the balance between cell and chemical potentials internally with the external chemical potential.As expected, given the same equilibrium conditions exist for both 1D and thin film models, the simulations reach the same equilibrium values for θ * and n * given the same parameter values and initial conditions.Interestingly, in deriving the thin film ODE model, we find that the scaled thin film length and height are equal.This relationship was assumed in an ad hoc manner in the work of Stevenson et al. (2010) when calculating cell traction.Our model thus provides a validation for their assumption in the case of uniform initial conditions.
For non-uniform initial conditions, we have derived small-time solutions, allowing for predictions of the stability of steady states.We also performed numerical simulations which show more complex behaviour compared to the spatially uniform case.The behaviour of the gel is still primarily driven by flow in the thin direction; however, the terms involving drag and spatial derivatives in the dependent variables are now active.Accordingly, we see spatial variations arising in θ p , n and h as the gel evolves.We have found that solutions exist where the gel height is non-uniform at steady state, even though the polymer fraction and cell density are spatially constant.When the initial polymer fraction θ I is non-uniform, local variations in the polymer induce solvent flow into or out of the gel to even out the polymer fraction; this correlates with spatial variations in the gel height at these points, reflecting the variations in mass.Similarly, when the initial gel height, h I , is non-uniform, then h will be spatially-varying to ensure θ p remains constant.
A significant difference in the thin geometry from the 1D case is that increasing drag now reduces spatial changes in the velocity.This is evident from equation (6.1d),where we see that increasing ξ decreases the influence of the spatial derivative terms on the velocity.A consequence of this is evident in Fig. 7a, which shows that, with large drag, initial spatial variations persist longer through the gel's evolution, as it is easier for fluid to flow vertically out of the gel than across the spatial domain due to the shearing forces present with a large drag coefficient.In the 1D case, we saw that increasing the drag coefficient tended to induce greater spatial gradients in the polymer and cells, as solvent could only enter and leave the gel through the endpoint at X = 1, and so had to flow across the entire spatial domain.We believe that this may explain why, for the thin film geometry, we have not observed oscillation between swelling and contraction, which we previously observed in the 1D geometry (Reoch et al. 2022).
We have noted that our thin film model derivation assumes that the non-dimensional diffusion coefficient D is O(1).Under this assumption, as seen in Sect.4, we cannot find equilibria in the thin film system with non-uniform polymer or cell distributions.Diffusion causes the cells to spread until a uniform cell density is reached across the gel.We note that examples with diffusion D = 0 in the 1D case also resulted in uniform equilibria.We have presented simulations with zero diffusion in the thin film here, finding examples where non-uniform equilibria persist in θ p and n.While we cannot definitively say that n = n(x, t) for the zero diffusion case, we have shown that, if such y-independent solutions can exist for D = 0, then we can find non-uniform equilibria in the thin film environment as well.We note that with small diffusion (D ≈ 0.0001), we do see quasi-steady states in θ p and n, where non-uniform states are found which slowly drift towards uniformity as a result of diffusive flux.
Unlike the studies by Trinschek et al. (2016Trinschek et al. ( , 2017) ) discussed in Sect. 1, we do not see situations where, at a certain point, the gel continues expanding lengthwise while its vertical swelling has stopped.Indeed, throughout these simulations, we see that the scaled height (or its mean value) is equal to the scaled gel length at equilibrium.As noted above, this supports modelling assumptions used in Stevenson et al. (2010), although we do note that, in the spatially non-uniform case, the average gel height h does not equal L for all time.In our model, there is no mechanism to cause a change in behaviour and decouple the lengthwise expansion from the vertical movement of the gel.This could possibly occur if the domain was somehow vertically or horizontally limited, or if an additional external pressure was imposed from one particular direction.A different cell force function might also see more emphasis on horizontal forces.An avenue for future work may be to explore whether these dynamics can be factored in, allowing for the expansion or contraction of a gel to be tailored in a certain direction, not just a situation where the length and height evolve similarly as seen here.
There is significant scope to extend this research and validate the model's behaviours through experimental collaboration.This would allow for parameter values to be fitted and suggest particular regions in the parameter space for deeper analysis to be carried out.The numerical results in this paper suggest further avenues to investigate experimentally, for example, confirming whether gels which are initially uniform in space retain this uniformity as they evolve, and whether small variations in the initial polymer profile do indeed result in spatially varying height profiles.
A number of modelling extensions are also of interest.Inclusion of surface tension would enable a study of its effect in cases where the gel height is found to vary; the aim would be to establish whether surface tension is sufficiently large to smooth out the height in such cases.Typical experimental setups, suggest future consideration of a thin gel on a solid substrate; this removes the symmetry about y = 0 so making the modelling more complex.Further, there may be cases where mechanical changes in the gel occur on a timescale similar to that of cell proliferation and death, making this an interesting scenario to explore.Again, this would add complexity to the model and its solution, requiring addition of suitable source terms to the PDE for cell density, and we have left this to future work.

B Thin film approximation: the y-independence of the first-order correction terms
In this appendix, we demonstrate that the O(ε) terms in cell density n 1 , pressure P 1 , polymer fraction θ p 1 , and polymer axial velocity v p 1 are also independent of y.This simplifies the derivation of an equation for the leading-order velocity (see Sect. 3.3).The result follows from essentially the same arguments as used in Sect.3.1.We begin by integrating equation (2.19) at O(ε −1 ) and using the no-flux boundary condition ∂n 1 /∂ y = 0 at y = h 0 , we find that n 1 = n 1 (x, t).This indicates that the first order correction to the cell force function is also independent of y, that is, G 1 = G 1 (x, t).Evaluating (2.22b) at O(ε −1 ) yields where 1 = P 1 − θ p 0 G 1 − θ p 1 G 0 .Using (2.25), we find that at y = h 0 , and hence 1 = 0. Therefore, Thus, the solvent axial velocity must also be independent of y at first order, i.e. v s 1 = v s 1 (x, t).