Mechanical cell interactions on curved interfaces

Pascal R Buenzli^a,[1] , Shahak Kuba^a, Ryan J Murphy^b, Matthew J Simpson^a

We propose a simple mathematical model to describe the mechanical relaxation of cells within a curved epithelial tissue layer represented by an arbitrary curve in two-dimensional space. The model represents the mechanics of the cell body either by straight springs between points of the curve, or by curved springs whose shape follows the curve. To understand the collective behaviour of these discrete models of cells at the broader tissue scale, we devise an appropriate continuum limit in which the number of cells is constant but the number of springs tends to infinity. The continuum limit shows that (i) the straight spring model and the curved spring model converge to the same dynamics; and (ii) the density of cells becomes governed by a diffusion equation in arc length space with second-order accuracy, where diffusion may be linear or nonlinear depending on the choice of the spring restoring force law. Our derivation of the continuum limit justifies that to reach consistent dynamics as the number of springs increases, the spring restoring force laws must be rescaled appropriately. Despite mechanical relaxation occurring within a curved tissue layer, we find that the curvature of the tissue does not affect tangential stress nor the mechanics-induced redistribution of cells within the layer in the continuum limit. However, the cell's normal stress does depend on curvature due to surface tension induced by the tangential forces. By characterising the full stress state of a cell, this model provide a basis to represent further mechanobiological processes.

[1]Corresponding author. Email address:

§ INTRODUCTION

Epithelial tissues are composed of confluent arrangements of cells that interact mechanically with each other through contact-induced cell body deformations [Ladoux & Mège, 2017]. Linking the mechanical properties of a tissue to those of its cells can help us understand how tissues grow and take their shape under spatial constraints as they interact with other tissues. It can also help us infer cell properties from experimentally observed tissue behaviour. Cells respond to mechanical cues, which may influence their proliferative behaviour, differentiation and survival [Opas, 1989, Weinans & Prendergast, 1996, Nelson et al., 2005, Keefer & Desai, 2011, Moeendarbary & Harris, 2014, Xi et al., 2019, Nelson, 2022]. The development of detailed mathematical models of cell and tissue mechanics is important to understand how cells behave in response to mechanics in morphogenesis and developmental biology, tumour growth, wound healing, and tissue engineering [Roose et al., 2003, Poujade et al., 2007, Basan et al., 2010, Murray et al., 2010, Tse et al., 2012, Huang,etal, 2014, Xi et al., 2019, Flegg et al., 2020, Murphy et al., 2021, Karakaya et al., 2022, Jensen & Revell, 2023].

Tissue mechanics is intricately related to tissue geometry. The geometric arrangement of cells within a tissue affects their mechanical interactions. Numerous studies have shown that geometry has a strong influence on tissue growth rate both in-vivo and in tissue engineering constructs [Rumpler et al., 2008, Bidan et al., 2012, Bidan et al., 2012, Bidan et al., 2016, Alias & Buenzli, 2017, Alias & Buenzli, 2018, Callens et al., 2020, Callens et al., 2023, Fratzl et al., 2022, Schamberger et al., 2023]. This geometric control of tissue growth is in part attributed to geometry-induced changes in cell mechanical stress which in turn influence cell behaviours through mechanobiological processes [Nelson et al., 2005, Luciano et al., 2021, Zmurchok, Bhaskar, & Edelstein-Keshet, 2018, Tambyah et al., 2020, Murphy et al., 2021]. Generally, it can be expected that where and how fast biological tissues grow and develop, is a combination of spatial constraints due to limited space availability, mechanical interactions, as well as biological cell behaviour [Dunlop et al., 2010, Gamsjäger et al., 2013, Goriely, 2017, Wang et al., 2018, Ambrosi et al., 2019].

In the present work, we focus on developing a mathematical model of mechanical interactions within a constant population of cells arranged in a confluent layer on a curved interface in two-dimensional space. We do not take into account other cell behaviours than the mechanical interactions of the cells. We aim to characterise how the distribution of cells along the curved interface evolves in space and time based on specifying the cells' mechanical properties, such as their stiffness and resting length. A key goal of this study is to understand whether the curvature of the interface influences this evolution as well as cellular stress, that can arise for example due to surface tension of the cellular tissue.

Our approach first considers two variations of a discrete, individual-based model of the mechanical behaviour of the cells. We then find how the evolution of cell density in these models can be captured, in an appropriate continuum limit, by a partial differential equation for a continuous cell density field. Both the consideration of discrete models and of their continuum limit is of strong interest in mathematical biology to bridge the cellular scale and the tissue scale. Discrete models allow us to represent the cellular entities that compose biological tissues and to encode cell behaviours directly. Continuum models capture average behaviour at the tissue scale more readily. However, without appropriate mathematical theories linking the cellular scale to the tissue scale, it can remain challenging to understand how cellular scale parameters relate to tissue scale behaviours, and to quantify these behaviours [Chopard & Droz, 1998, Li et al., 2022, VandenHeuvel, Buenzli & Simpson, 2024].

Many models of cellular dispersion are based on random walks [Codling, Plank & Benhamou, 2008]. These types of models are appropriate when motile cells migrate through free space and their interactions represent crowding constraints [Chopard & Droz, 1998, Simpson et al., 2010, Li, Buenzli & Simpson, 2022, Li et al., 2022], but they are not well-suited to modelling mechanical interactions occurring within confluent layers of cells. Previous works have studied mechanical interactions between confluent cells arranged along a one-dimensional axis [Murray et al., 2009, Murray et al., 2011, Murray et al., 2012, Fozard et al., 2010, Murphy et al., 2019, Murphy et al., 2020, Murphy et al., 2021, Tambyah et al., 2020]. These works find that when cells interact mechanically like overdamped Hookean springs with stiffness constant $k$ and drag coefficient $\eta$, the density of cells $q(x,t)$ at position $x$ along the line and time $t$ evolves, in an appropriate continuum limit, according to the nonlinear diffusion equation

\begin{align}\label{cell-density-1D}
\pd{q}{t} = \pd{}{x}\left(\frac{k}{\eta}\frac{1}{q^2}\pd{q}{x}\right).
\end{align}

This partial differential equation represents the mechanical relaxation of the tissue as a fast diffusion process with diffusivity $D(q) = k/(\eta q^2)$, in which low densities disperse rapidly [Newman, 1980, King & McCabe, 2003]. It belongs to a class of porous medium equations with a negative exponent $-2$ for the density dependence of the diffusivity, representing fast diffusion [Vázquez, 2006].

Equation (<ref>) provides a direct relationship between mechanical properties of a cell body ($k$ and $\eta$) and the mechanics-driven evolution of a continuous cell density at the tissue scale. For non-Hookean restoring force laws, other nonlinear density dependences of the diffusivity hold [Murray et al., 2012, Baker, Parker & Simpson, 2019]. In this paper, we generalise Eq. (<ref>) to represent the evolution of cell density due to mechanical relaxation in curved epithelial layers. We base our mathematical model on a discrete model of cells interacting as a chain of mechanical springs similar to those proposed in the above works, except that we consider the springs to be oriented along a curved interface. This interface represents a substrate providing support for the cells, such as an epithelial basement membrane, or other biological tissues supporting confluent cell layers, such as osteoblasts lying on bone tissue [Martin, Burr & Sharkey, 1998]. We propose two different spring models to represent the cells' mechanical interactions: one in which the springs are straight and bridge across discrete positions of the continuous interface, and one in which the springs are curved and follow the curved interface. We then devise a new continuum limit procedure of these spring models of cellular mechanics based on assuming that each cell comprises $m$ springs and letting $m\to\infty$. In this way, the number of cells and cell density can remain finite while the length of the springs goes to zero, allowing us to perform expansions in terms of small spring lengths. Our procedure allows us to provide formal justifications for cell parameter rescalings with $m$ that are necessary for converging to consistent dynamics. These scalings were previously mentioned in [Murray et al., 2009, Murphy et al., 2019] without rigorous justification. Our procedure also provides the order of approximation of the continuum limit to describe discrete model behaviour as $m\to\infty$, which explains why continuous densities approximate cell densities from the discrete model particularly well at the midpoint of the springs. We also investigate time scales of mechanical relaxation in the discrete models and the stability of homogenous densities depending on the spring restoring force laws [Murray et al., 2012]. Finally, we calculate both the tangential and normal stresses within the cellular tissue to illustrate how surface tension induced by the spring forces generates curvature-dependent normal stress.

§ DISCRETE MODEL OF MECHANICAL CELL INTERACTIONS

Following previous works [Murray et al., 2009, Murray et al., 2011, Murray et al., 2012, Fozard et al., 2010, Baker, Parker & Simpson, 2019, Murphy et al., 2019, Murphy et al., 2020, Murphy et al., 2021, Tambyah et al., 2020], we consider mechanical interactions between adjacent cells, which we interpret as a model of the cross-section of an epithelial cell monolayer (Figure <ref>). Each cell is in contact with a single cell on either side. While these previous works assumed cells to be constrained to the $x$-axis, here we consider that cells lie on a static interface in two-dimensional space. The interface is described by a general curve $\b r(s)=\big(x(s), y(s)\big)$ parametrised by arc length $s\in [0,L]$, where $L$ is the total length of the interface. The interface may be an open curve, or a closed loop with $\b r(0)=\b r(L)$. It can have arbitrary shape, except that we assume $\b r(s)$ to be differentiable, meaning that there are no sharp corners and a tangent vector exists everywhere. This assumption is physically realistic for epithelial layers, and mathematically convenient to describe the orientation of the cells and the direction of mechanical forces acting on them.

In the following, we first present two discrete models of mechanical interaction occurring between the cells on the interface, one based on straight springs, and one based on curved springs. We then show in Section <ref> that in an appropriate continuum limit in which the number of springs goes to infinity, both models converge to the same nonlinear diffusion equation for the local cell density.

Both discrete models assume that the shape and mechanical behaviour of a biological cell can be represented by a chain of $m$ internal mechanical springs along the cell body. These springs are connected to those of the neighbouring cells (Figure <ref>). To simplify notation, throughout the paper we denote by $i=0,\ldots, M$ the index of a node connecting two springs irrespective of whether this node is internal to a cell, or at the boundary between two cells. We denote the total number of cells by $N$, so that the total number of springs is $M=mN$. If the interface is a closed loop, the last spring is connected to the first spring, so that node $i=M$ is identified with node $i=0$ (there are only $M$ distinct nodes in this case).

Cell body motion occurs in dissipative environments where inertial effects can be neglected. We therefore assume an overdamped regime, in which the position $\b r_i(t)$ of node $i$ at time $t$ evolves according to

\begin{align}\label{discrete-evo-0}
\eta \td{\b r_i}{t} = \b F_i^{(-)} + \b F_i^{(+)} + \b F_i^{(\text{n})},
\end{align}

where $\eta$ is the viscous drag coefficient, and $\b F_i^{(-)}$ and $\b F_i^{(+)}$ are the mechanical forces exerted on node $i$ due to the springs connecting it to nodes $i-1$ and $i+1$, respectively. To ensure that the nodes remain on the interface at all times, an additional normal reaction force $\b F_i^{(\text{n})}$ may be required in curved geometries, which represents the interaction of the cells with the substrate (Figure <ref>). This additional force is not required in models of epithelia on flat interfaces [Murray et al., 2009, Murphy et al., 2019]. We now detail how the forces $\b F_i^{(-)}$, $\b F_i^{(+)}$, and $\b F_i^{(\text{n})}$ are defined in the straight spring model, and in the curved spring model.

(a) Straight spring model and (b) curved spring model of the mechanical interaction between cells along an interface $\b r(s)$ (solid black line). In this figure, each cell is composed of $m=2$ springs (gray coils). Nodes within a cell are shown as white circles, and nodes connecting two cells are shown as black circles. (a) The force diagram for the straight spring model shows that restoring forces are directed along the secant line between two nodes on the interface. The normal reaction force $\b F_i^{(\text{n})}$ ensures that the net force is parallel to the unit tangent vector $\b \tau_i$ at node $i$. (b) In the curved spring model, the restoring forces $\b F_i^{(\pm)}$ at node $i$ are already parallel to $\b \tau_i$ and so there is no normal reaction force $\b F_i^{(\text{n})}$.

§.§ Straight spring model

In the straight spring model, the line of action of the spring restoring force is directed along the secant line connecting two nodes (Figure <ref>a), i.e.,

\begin{align}
\b F_i^{(-)} = - f\big(\|\b r_i-\b r_{i- 1}\|\big) \frac{\b r_i-\b r_{i- 1}}{\|\b r_i-\b r_{i-1}\|}, \quad \b F_i^{(+)} = f\big(\|\b r_{i+1}-\b r_i\|\big) \frac{\b r_{i+1}-\b r_i}{\|\b r_{i+1}-\b r_i\|},
\end{align}

where $f(\ell)$ is the restoring force law that depends on spring length $\ell$. The normal reaction force $\b F_i^{(\text{n})}$ in Eq. (<ref>) is defined to be exactly opposite to the component of $\b F_i^{(-)}+\b F_i^{(+)}$ normal to the interface at node $i$, so that the net force and velocity $\d\b r_i/\d t$ of node $i$ in Eq. (<ref>) are tangent to the interface at all times. Denoting with $\b\tau_i$ the unit tangent vector to the interface at node $i$ in the direction of increasing index $i$, the position of node $i$ evolves in the straight spring model as

\begin{align}\label{straight-spring-evo-1}
\eta \td{\b r_i}{t} = \left[ \left( f\big(\|\b r_{i+1}-\b r_i\|\big)\frac{\b r_{i+1}-\b r_i}{\|\b r_{i+1}-\b r_i\|}  - f\big(\|\b r_i-\b r_{i-1}\|\big)\frac{\b r_i-\b r_{i-1}}{\|\b r_i-\b r_{i-1}\|}\right)\vdot \b\tau_i \right]\b \tau_i.
\end{align}

Since the nodes remain on the interface at all times, we can reference their position along the parametric curve $\b r(s)$ by their arc length coordinate $s_i(t)$, such that

\begin{align}\label{pos-vs-arclength}
\b r_i(t) = \b r\big(s_i(t)\big).
\end{align}

Using the fact that the unit tangent vector of the interface at $s$ is $\b\tau(s)=\d\b r(s)/\d s$ and that $\b\tau_i = \b\tau\big(s_i(t)\big)$, we obtain from Eq. (<ref>),
\begin{align}\label{velocity-vs-arclength}
\td{\b r_i(t)}{t} = \td{\b  r\big(s_i(t)\big)}{s}\td{s_i(t)}{t} = \b \tau_i\td{s_i(t)}{t}.
\end{align}

Substituting these expressions for $\b r_i$ and $\d\b r_i/\d t$ in terms $s_i$ and $\d s_i/\d t$ in Eq. (<ref>) provides the evolution of the nodes' arc length coordinates:

\begin{align}\label{straight-spring-evo-2}
\eta \td{s_i}{t} = \left[f(\|\b r(s_{i+1})-\b r(s_i)\|)\frac{\b r(s_{i+1})-\b r(s_i)}{\|\b r(s_{i+1})-\b r(s_i)\|} - f(\|\b r(s_i)-\b r(s_{i-1})\|)\frac{\b r(s_i)-\b r(s_{i-1})}{\|\b r(s_i)-\b r(s_{i-1})\|}\right]\vdot\b\tau_i.
\end{align}

Since $\b r(s)$ is a known parametric curve, we can solve Eq. (<ref>) instead of Eq. (<ref>) to evolve the discrete model, with appropriate boundary conditions when $i=0$ and $i=M$ (see Section <ref>).

§.§ Curved spring model

An alternative model of mechanical interactions is to consider curved springs that follow the curvature of the interface (Figure <ref>b), so that the length of a spring is the arc length between two nodes, and the line of action of the restoring forces $\b F_i^{(\pm)}$ is directed along the unit tangent $\b\tau_i$ at node $i$, i.e.,

\begin{align}
\b F_i^{(-)} = -f(s_i-s_{i-1})\b\tau_i, \qquad \b F_i^{(+)} = f(s_{i+1}-s_i)\b\tau_i,
\end{align}

where $s_i$ is the arc length coordinate of node $i$. Since these forces are tangent to the interface, $\b F_i^{(\text{n})}=\b 0$. In the curved spring model, the position of node $i$ therefore evolves from Eq. (<ref>) as

\begin{align}
\eta\td{\b r_i}{t} = \left[ f(s_{i+1}-s_i) - f(s_i-s_{i-1})\right]\b \tau_i.
\end{align}

Using Eq. (<ref>), the evolution of arc length coordinates in this model is given by

\begin{align}\label{curved-spring-evo-2}
\eta\td{s_i}{t} = f(s_{i+1}-s_i) - f(s_i - s_{i-1}).
\end{align}

Mathematically, Eq. (<ref>) is of the same form as the evolution of coordinates $x_i$ in a straight chain of springs along the $x$-axis [Murray et al., 2009, Murray et al., 2011, Murray et al., 2012, Murphy et al., 2019], except that $x_i$ are replaced by arc length coordinates $s_i$. It is interesting to note that Eq. (<ref>) is manifestly independent of the shape of the parametric curve $\b r(s)$, but Eq. (<ref>) is not. When nodes interact mechanically via straight springs, the evolution of their node positions along the interface depends on the shape of the interface.

§.§ Restoring force law

In both the straight spring model and the curved spring model we may consider a general restoring force law $f(\ell)$ that depends on the evolving length $\ell$ of the spring. In our numerical simulations, we consider three different force laws, including a Hookean restoring force law linear in the elongation of the spring, and two nonlinear restoring force laws. The Hookean restoring force law is

\begin{align}\label{Hookean-restoring-force}
f(\ell) = k(\ell-a),
\end{align}

where $a$ is the resting spring length and $k$ is the spring stiffness.

Hooke's law means that cells can in principle be squeezed to zero length since the force remains finite when $\ell=0$. It also means that the further apart cells are, the stronger their mechanical interaction. These properties are biologically unrealistic for very small $\ell$ and for very large $\ell$, so we also consider nonlinear restoring force laws that diverge when $\ell\to 0$ and that remain bounded when the separation between cells increases ($\ell\to\infty$). Our specific choices of such restoring force laws are such that in the continuum limit, the mechanical relaxation of cells will be described by linear diffusion of cell density for our first choice, and by a common type of porous medium diffusion for our second choice (see Section <ref>, Table <ref>). The first nonlinear restoring force law is

\begin{align}\label{nonlinear-restoring-force}
f(\ell) = ka^2\left(\frac{1}{a}-\frac{1}{\ell}\right),
\end{align}

where $\rho_0=1/a$ is the resting spring density. The scaling factor $ka^2$ of this restoring force is such that the linearisation of Eq. (<ref>) about $\ell=a$ is the Hookean restoring force (<ref>) with spring constant $k$ and resting length $a$ (Figure <ref>).

The second nonlinear restoring force law is

\begin{align}\label{porous-restoring-force}
f(\ell) = \frac{k a^3}{2}\left(\frac{1}{a^2} - \frac{1}{\ell^2}\right).
\end{align}

Likewise, the scaling factor $ka^3/2$ ensures that the linearisation of Eq. (<ref>) about the resting length $a$ matches the Hookean restoring force (Figure <ref>).

Comparison between the Hookean restoring force law $f(\ell)=k(\ell-a)$ (magenta), the nonlinear restoring force law $f(\ell)=ka^2(1/a-1/\ell)$ (green), and the nonlinear restoring force law $f(\ell)=(ka^3/2)\big(1/a^2-1/\ell^2\big)$ ($a=1$, $k=1$ in arbitrary units). The nonlinear restoring force scaling factors are such that their linearisation about the resting spring length $\ell=a$ gives the Hookean restoring force.

The restoring force laws $f(\ell)$ in Eqs (<ref>)–(<ref>) are defined to be positive when the spring is elongated, and negative when the spring is compressed. We can define tangential stress along the cell body based on these force laws as follows. Considering first Hooke's law, the spring constant $k$ of a portion of cell body of resting length $a$ and cross-sectional area $A$ is $k=EA/a$, where $E$ is the cell body's Young modulus. The tangential stress along the cell body can thus be defined at each inner spring of a cell by

\begin{align}\label{tangential-stress}
\sigma_{\tau\tau}(\ell) = - \frac{f(\ell)}{A} = - \frac{E}{k}\frac{f(\ell)}{a},
\end{align}

with the sign convention that compressive stress is positive and tensile stress is negative. For all the restoring forces that we consider, we will calculate tangential stress as $\sigma_{\tau\tau}/E$ using Eq. (<ref>) as it only depends on $k$, $a$, and $\ell$ and it is always proportional to stress. The quantity $\sigma_{\tau\tau}/E$ represents strain for Hookean springs but not for nonlinear springs.

§.§ Numerical simulations

Numerical simulations of the discrete model are based on evolving arc length coordinates using Eq. (<ref>) for straight springs and Eq. (<ref>) for curved springs, supplemented with the following boundary conditions. If the interface is an open curve, we fix the position of the boundary nodes $i=0$ and $i=M$ at the curve's end points, so that their arc length coordinates satisfy

\begin{align}\label{fixed-bc}
s_0(t) = 0, \qquad s_M(t) = L, \qquad t\geq 0.
\end{align}

If the interface is a closed loop, we consider periodic boundary conditions on the node indices by identifying node $i=M$ with node $i=0$, so that

\begin{align}\label{pbc}
s_0(t)=s_M(t), \qquad t\geq 0.
\end{align}

Our numerical simulations consider several interfaces for which we require an arc length parametrisation. Finding arc length parametrisations usually involves reparametrising an interface described by another, known parametrisation. For a circular interface of radius $R$ parametrised by an angle $\theta\in[0,2\pi)$, arc length is $s=R\theta$ so $\b r(s) = \big(R\cos(s/R), R\sin(s/R)\big)$, and the unit tangent vector is $\b\tau(s) = (-\sin(s/R), \cos(s/R))$. For other interface shapes, arc length $s$ is calculated from a known parametrisation $\widetilde{\b r}(u)$ of the interface by numerically integrating

\begin{align}
s(u) = \int_0^u\!\!\!\d u\ \|\widetilde{\b r}'(u)\| \label{arclength}
\end{align}

using the composite Simpson's 1/3 rule. The arc length parametrisation of the interface is then obtained as $\b r(s)=\widetilde{\b r}\big(u(s)\big)$ where $u(s)$ is the inverse function of $s(u)$. This inverse function is estimated numerically by finding $u$ in Eq. (<ref>) such that $s(u)-s=0$ using the root-finding algorithm  of 's standard library, which is based on the  algorithm [LDC2 compiler, 2024]. All our simulations assume arbitrary units of space and time. These units can be re-scaled to match any particular application by working with a relevant length scale defined by the cell diameter, and a mechanical relaxation time scale determined by the ratio $\eta/k$ (Section <ref>).

Equations (<ref>) and (<ref>) are solved numerically using a simple explicit forward Euler time-stepping scheme with time step $\Deltaup t=0.001$, and with boundary conditions (<ref>) for open curves, and (<ref>) for closed curves. Initial conditions assume that all cell boundaries are evenly spaced except for one that is offset by half a resting cell length. Spring boundaries are then initialised to be evenly spaced within a cell. For full detail on the algorithms used to solve the discrete models, the reader is referred to the  computer code available on  [Buenzli et al., 2024].

Time snapshots of the mechanical relaxation of $N=4$ cells with $m=4$ inner springs (stress-coloured coils) along the open curve $\widetilde{\b r}(u)=\big(u,R\sin(u)\big)$ (solid black curve) using the curved spring model and a Hookean restoring force. Cell boundaries are shown as black circles. Inner spring boundaries are shown as open circles. The resting spring length is chosen such that the steady state is stress-free; $R=0.8$, $k=4$, $\eta=0.25$, $a\approx 0.45$, $\Deltaup t = 0.001$.

We start by considering a sinusoidal interface given by $\widetilde{\b r}(u)=\big(u, R\sin(u)\big)$ for $u\in[0,2\pi]$ (Figure <ref>). The interface is populated with $N=4$ cells, each containing $m=4$ curved, Hookean springs, so that there are $M=16$ springs in total. The first and last spring boundaries are fixed at $\widetilde{\b r}(0)=(0,0)$ and $\widetilde{\b r}(2\pi)=(2\pi,0)$, respectively. The spring resting length is chosen as $a=L/M\approx 0.45$, where $L\approx 7.2$ is the length of the curve between $u=0$ and $u=2\pi$. Thus, springs are stress-free when in mechanical equilibrium. Since $m=4$, the resting cell length is $4a\approx 1.8$. In the initial configuration, the boundary between the second and third cell is displaced, such that the second cell is elongated by half a resting cell length, and the third cell is compressed by a half a resting cell length. Time snapshots of the numerical simulation are shown in Figure <ref> where each spring $i$ of length $\ell_i=(s_i-s_{i-1})$ is coloured according to its normalised tangential stress $\sigma_{\tau\tau}(\ell_i)/E$. These snapshots show that the springs relax mechanically to even out mechanical stress by redistributing the position of their extremities along the curved interface.

The evolution of the arc length positions of the spring nodes for the same simulation is shown in Figure <ref> (solid black and grey lines). In Figure <ref>a, each spring $i$ is coloured according to its normalised tangential stress $\sigma_{\tau\tau}(\ell_i)/E$. In Figure <ref>b, each spring $i$ is coloured according to its local cell density $q_i=1/(m\ell_i)$. By evening out tangential stress along the interface with time, the mechanical forces between the cells also tend to even out cell density along the interface.

Evolution of spring boundary positions (thin gray lines) and cell boundary positions (thick black lines) in the simulations shown in Figure <ref>. Each spring is coloured according to (a) its tangential stress $\sigma_{\tau\tau}(\ell_i)/E$; and (b) cell density $q_i$.

§ CONTINUUM LIMIT

In both the straight spring model and curved spring model, the cells relax to an equilibrium state through their mechanical interactions. The collective relaxation of the cells corresponds to the mechanical relaxation of the tissue. This tissue-scale relaxation can be described in a continuum limit by deriving the conservation equation that governs the evolution of cell density when the number of cells $N$ is kept constant and the number of springs per cell $m$ goes to infinity, $m\to\infty$ [Murphy et al., 2019]. This way of taking the continuum limit allows cell density to remain finite on a finite interface length $L$, while the total number of springs in the system goes to infinity, $M\to\infty$.

We first show in Section <ref> that as $m\to\infty$ the dynamics of the straight spring model converges to that of the curved spring model. In Section <ref>, we then propose an alternative and more intuitive derivation of the continuum limit of the evolution of cell density than the derivations provided in Murray et al., 2009, Murphy et al., 2019, Tambyah et al., 2020. This new derivation enables us to find how the drag coefficient $\eta$, the spring stiffness $k$, and the resting length $a$ must be rescaled for the dynamics of the system to converge as the number of springs per cell increases. Furthermore, in Section <ref> we propose an alternative justification for these parameter scalings based on analysing time scales of mechanical relaxation in the discrete models.

§.§ Equivalence of straight and curved springs in the continuum limit

When the number of springs per cell $m$ goes to infinity, the distance between the nodes tends to zero, so that $\Deltaup s = s_{i+1}-s_{i}\to 0$. Using Eq. (<ref>) and the differentiability of $\b r(s)$, we therefore have

\begin{align}
&\b r_{i+1} = \b r(s_i + \Deltaup s) = \b r_i + \b\tau_i\Deltaup s + \order(\Deltaup s), \notag
\\&\b r_{i-1} = \b r(s_i-\Deltaup s) = \b r_i-\b\tau_i\Deltaup s + \order(\Deltaup s), \label{contlim-pos}
\end{align}

where $\order(\Deltaup s)$ terms are such that $\lim_{\Deltaup s\to 0}\order(\Deltaup s)/\Deltaup s=0$ [Olver, 1974]. Equations (<ref>) show that the Euclidean distance between two nodes converges to the arc length distance along $\b r(s)$,

\begin{align}\label{contlim-arclength}
\|\b r_{i+1}-\b r_i\| \sim \Deltaup s = s_{i+1}-s_i, \qquad m\to \infty,
\end{align}

and that the lines of force of $\b F_i^{(\pm)}$ in the straight spring model become parallel to the unit tangent vector $\b \tau_i$ at $\b r_i$:

\begin{align}\label{contlim-tangent}
&\frac{\b r_{i+1}-\b r_{i}}{\|\b r_{i+1}-\b r_{i}\|} \to \b\tau_i, \qquad\frac{\b r_{i}-\b r_{i-1}}{\|\b r_{i}-\b r_{i-1}\|} \to \b\tau_i, \qquad m\to\infty.
\end{align}

Using Eqs (<ref>)–(<ref>), the evolution equation (<ref>) for arc length coordinates in the straight spring model converges to the evolution equation (<ref>) of the curved spring model as $m\to\infty$:

\begin{align}
\eta\td{s_i}{t} &\sim \big(f(s_{i+1}-s_i)\b\tau_{i} - f(s_i-s_{i-1})\b\tau_i\big)\vdot\b\tau_i,\notag
\\&= f(s_{i+1}-s_i)-f(s_i-s_{i-1}), \qquad m\to\infty. \label{discrete-evo}
\end{align}

Thus, for large number of springs per cell $m$, the evolution of spring boundaries for straight springs tends to that of curved springs.

Evolution of spring boundary positions along the circle $\b r(s)=\big(R\cos(s/R),R\sin(s/R)\big)$ with straight and curved springs, $N=4$, $m=1$, $R=1$, $k=1$, $\eta=1$, $\Deltaup t=0.001$. (a) Initial condition with straight springs of resting length $a=2\sin(\pi/4) = \sqrt{2}$; (b) Initial condition with curved springs of resting length $a=2\pi/4$; (c) Comparison of the dynamics of mechanical relaxation between straight and curved springs.

To illustrate the difference between the straight spring model and the curved spring model, we present numerical simulations in Figure <ref> with few springs, $N=4$, $m=1$, and a circular interface. The spring resting length is chosen such that springs are stress-free when in mechanical equilibrium, meaning that $a=2\pi R/M$ for curved springs, and $a=2R\sin\big(\pi/M\big)$ for straight springs, which is the edge length of an $M$-sided regular polygon inscribed in a circle of radius $R$. When one of the nodes is offset by an arc length distance $a/2$ along the circle, resulting in one compressed cell adjacent to one extended cell, the dynamics of mechanical relaxation is slightly different between the straight spring model and the curved spring model (Figure <ref>c). However, the mechanical equilibrium state is identical.

Figure <ref> shows that if the interface is more complex, both the dynamics and the equilibrium state between straight springs and curved springs might not be the same. In the straight spring model, some nodes can remain stuck in regions of the interface with high curvature, due to the particular arrangement of spring force directions (Figure <ref>a). In contrast, spring force directions are always tangent to the interface with curved springs, so nodes can always slide away from regions with high curvature (Figure <ref>b). These differences between the straight spring and curved spring models quickly disappear as $m$ is increased. Figure <ref> compares the arc length positions of the spring nodes when $m=1$ and when $m=8$ in both spring models. It shows that both the dynamics and steady state of mechanical relaxation on the cross-shaped boundary of Figure <ref> become indistinguishable already for $m=8$ springs per cell (Figure <ref>b). However, to obtain similar dynamics and relaxation times of the spring nodes when $m=1$ and when $m=8$, it is necessary to rescale the spring resting length by a factor $1/m$, and the ratio $k/\eta$ by a factor $m^2$, which we do by rescaling $\eta$ by $1/m$ and $k$ by $m$. While these scalings were already proposed by Murray et al., 2009, Murphy et al., 2019, there were not previously justified mathematically. We will justify these scalings rigorously in Sections <ref> and <ref>.

Comparison of mechanical relaxation between straight and curved springs on a cross-shaped interface with $N=8$ cells and $m=1$ spring per cell. The interface is defined in polar coordinates by the polar equation $R(\theta)=R_0\big(\cos^4(\theta) + \sin^4(\theta)\big)$;  $R_0=1$, $\eta=1$, $k=1$, $\Deltaup t=0.001$. The resting length $a$ is chosen such that there is no tangential stress in steady state. The snapshots show the initial configuration ($t=0$) and a mechanically relaxed configuration ($t=100$) for (a) straight springs with $a\approx 0.7368$; and (b) curved springs with $a\approx 0.8010$.

Comparison of the evolution of spring boundary positions along the cross-shaped interface of Figure <ref> between straight springs and curved springs with $N=8$ cells. Cell boundaries are shown as solid black line (straight spring model) and thick dashed green lines (curved spring model). Inner spring boundaries within the cells are shown as thin grey lines (straight spring model) and thin dashed green lines (curved spring model). Spring resting lengths are chosen such that there is no tangential stress in steady state. (a) $m=1$ spring per cell, $k=1$, $\eta=1$, $a \approx 0.7368$ for straight springs, $a\approx 0.8010$ for curved springs, $\Deltaup t = 0.001$; (b) $m=8$ springs per cell, $k=8$, $\eta=1/8$, $a\approx 0.0995$ for straight springs, $a\approx 0.1001$ for curved springs, $\Deltaup t = 0.001$.

§.§ Evolution of cell density in the continuum limit

Equation (<ref>) is identical to the evolution of node coordinates $x_i$ in a straight chain of springs along the $x$ axis [Murray et al., 2009, Murphy et al., 2019], with the only difference that the coordinates are now arc lengths $s_i$ along a curved interface. The geometry of the interface does not appear explicitly in this equation. The evolution of cell density in the continuum limit is therefore identical to that found in [Murray et al., 2009, Murphy et al., 2019] subject to substituting $x$ for arc length $s$. Murray et al. [Murray et al., 2009, Murray et al., 2011, Murray et al., 2012] and Tambyah et al. [Tambyah et al., 2020] derive the continuum limit based on considering density as the reciprocal of the metric $\p x_i/\p i$ for a continuous spring index $i$. Murphy et al. [Murphy et al., 2019, Murphy et al., 2020] derive the same continuum limit based on local spatial averages. We provide an alternative derivation, based on a conceptually simpler expansion of cell density in small spring lengths. This new derivation has the benefit of (i) clearly justifying scalings required on parameters of the discrete model to reach consistent dynamics as the number of springs increases; (ii) calculating the order of approximation of the discrete model provided by the continuum limit. This order of approximation justifies in particular why in numerical simulations the continuous density matches the discrete density near the midpoint between spring nodes, which was previously observed, but not explained [Murphy et al., 2019].

Since the straight spring model becomes identical to the curved spring model when $m\to\infty$, we consider the curved spring model and denote by $\ell_i = (s_i-s_{i-1})$ the length of the $i^\text{th}$ spring. The local spring density in the discrete model is given by $\rho_i = 1/\ell_i$.
The evolution of spring density is calculated from Eq. (<ref>) as:

\begin{align}
\td{\rho_i}{t} &= -\frac{1}{\ell_i^2}\td{\ell_i}{t} = -\frac{1}{\ell_i^2}\left(\td{s_i}{t}-\td{s_{i-1}}{t}\right)\notag
\\&= -\frac{1}{\eta\ell_i^2}\Big( f(s_{i+1}-s_i) - f(s_i-s_{i-1}) - \big(f(s_i-s_{i-1}) - f(s_{i-1}-s_{i-2})\big)\Big)\notag
\\&= -\frac{1}{\eta}\frac{f(\ell_{i+1}) - 2f(\ell_i) + f(\ell_{i-1})}{\ell_i^2}.\label{evo-spring-density-discrete}
\end{align}

Equation (<ref>) is exact for curved springs, and asymptotic for straight springs when $m\to\infty$ via Eq. (<ref>). To take the continuum limit, we now introduce a continuous function $\rho(s,t)$ of arc length position $s$ and time $t$, such that

\begin{align}\label{rho-cont}
\rho\big(\overline{s}_i, t\big) = \rho_i(t) = 1/\ell_i(t),
\end{align}

\begin{align}\label{midpoint-arclength}
\overline{s}_i  = \frac{1}{2}\big(s_i + s_{i-1}\big)
\end{align}

is the (time-dependent) arc length midpoint between the spring boundaries $s_{i-1}$ and $s_i$. We also define the restoring force law as a function of spring density instead of length by introducing

\begin{align}\label{f-tilde}
\widetilde{f}(\rho) = f(1/\rho).
\end{align}

Because spring lengths tend to zero, spring density diverges as $\Order(m)$ as $m\to\infty$. Since each cell contains $m$ springs, we also introduce a continuous, local cell density function $q(s,t)$ of arc length position $s$ and time $t$ based on $\rho(s,t)$, defined as

\begin{align}\label{q}
q(s,t) = \frac{\rho(s,t)}{m}.
\end{align}

We thus expect cell density $q(s,t)$ to remain finite, of order $\Order(1)$, in the continuum limit. However, we will see that additional requirements on the restoring force and drag coefficient are necessary for the evolution of $q(s,t)$ in the continuum limit to match the evolution of the discrete models with $m \gg 1$. We first explore the large-$m$ behaviour of Eq. (<ref>) for the (diverging) spring density $\rho(s,t)$, and then consider under what conditions the evolution of cell density $q(s,t)$ is well defined in the limit $m\to\infty$.

Using the continuous function $\rho(s,t)$, Eq. (<ref>), and the fact that $\ell_i = (s_{i}-s_{i-1}) = 1/\rho(\overline{s}_i,t)$, Eq. (<ref>) can be rewritten as

\begin{align}\label{evo-density-midpoint}
\td{}{t}\rho(\overline{s}_i,t)  &\sim -\frac{1}{\eta}\frac{\widetilde{f}\big(\rho(\overline{s}_{i+1},t)\big)-2\widetilde{f}\big(\rho(\overline{s}_i,t)\big) + \widetilde{f}\big(\rho(\overline{s}_{i-1},t)\big)}{\ell_i^2}, \qquad m\to\infty.
\end{align}

The right hand side is similar to a discretised second order derivative of $\widetilde{f}\big(\rho(s,t)\big)$ with respect to arc length evaluated at $\overline{s}_i$. However, $\overline{s}_{i-1}$, $\overline{s}_i$, and $\overline{s}_{i+1}$ are not evenly spaced, and the time derivative in the left hand side must also take into account that $\overline{s}_i$ is time dependent. In Appendix <ref>, we show that by properly expanding $\overline{s}_{i+1}$ and $\overline{s}_{i-1}$ about $\overline{s}_i$ in the limit $m\to\infty$ and accounting for the time dependence of $\overline{s}_i$ in the left hand side, Eq. (<ref>) asymptotically becomes the partial differential equation

\begin{align}
\pd{\rho}{t}(\overline{s}_i,t) &= -\frac{1}{\eta}\pd{^2}{s^2}\widetilde{f}\big(\rho(\overline{s}_i,t)\big) \left(1+ \Order\left(\frac{1}{m^2}\right)\right),\quad m\to\infty.\label{evo-spring-density-continuous}
\end{align}

It is important to emphasise that the $\Order(1/m^2)$ correction relies on matching the continuum and discrete spring densities at the arc length midpoint of a spring in Eq. (<ref>). Matching at other points along the spring decreases accuracy to $\Order(1/m)$ (Appendix <ref>). As expected from these results, the arc length midpoint of a spring is precisely where numerical simulations show excellent agreement between discrete densities and continuum densities (Section <ref>). The order of accuracy of the continuum limit was not known in previous works.

The evolution of cell density can now be obtained from Eq. (<ref>) by substituting $\rho(s,t)=m q(s,t)$. This gives, for an arbitrary arc length coordinate $s$:

\begin{align}
\pd{q}{t}(s, t) = -\pd{^2}{s^2}\left(\frac{1}{\eta m}\widetilde{f}\big(m q(s,t)\big)\right) + \Order\left(\frac{1}{m^2}\right) \label{evo-cell-density-1}.
\end{align}

For this partial differential equation to be well defined as $m\to\infty$, the restoring force and drag coefficient must be such that the following limit exists:

\begin{align}\label{restoring-force-contlim}
\lim_{m\to\infty}\frac{1}{\eta m}\widetilde{f}\big(m q(s,t)\big) = F\big(q(s,t)\big),
\end{align}

where the limit defines a cell-density-dependent restoring force function $F(q)$. If this limiting force function $F(q)$ is reached with order $\Order(1/m^2)$, then the evolution equation for cell density in the continuum limit is a conservation equation with flux given by the gradient of $F$ with respect to arc length $s$, $\p F\big(q(s,t))/\p s$, up to second-order corrections in $1/m$:

\begin{align}\label{evo-cell-density-2}
\pd{q}{t}(s,t) = -\pd{^2}{s^2}F\big(q(s,t)\big) + \Order
\left(\frac{1}{m^2}\right).
\end{align}

We can recast the limiting equation as the nonlinear diffusion equation

\begin{align}\label{nonlinear-diffusion}
\pd{q}{t}(s,t) = \pd{}{s}\left(D\big(q(s,t)\big)\pd{q}{s}(s,t)\right), \qquad \text{where}\quad D(q) = - \td{F}{q}(q).
\end{align}

Equation (<ref>) was derived previously for arbitrary restoring force laws by Murray et al., 2012 but for the diverging density of springs, and therefore without providing rigorous conditions on the required scalings of the restoring force law with $m$ exhibited in Eq. (<ref>). Equation (<ref>) was derived by Murphy et al., 2019 for the cell density, which is well-defined in the continuum limit, when the restoring force is Hookean, but without providing rigorous justifications for the restoring force law scaling and the order of approximation.

Solutions to the partial differential equations (<ref>) or (<ref>) require the specification of initial and boundary conditions, which are chosen to match those of the discrete model. When the interface is an open curve, fixing the position of node $i=0$ and $i=M$ in the discrete model in Eqs (<ref>) means that no spring can move past the fixed boundary nodes $\b r_0(t) = \b r(0)$ and $\b r_M(t) = \b r(L)$. Accordingly, no-flux boundary conditions on $q(s,t)$ are imposed at $s=0$ and $s=L$:

\begin{align}\label{no-flux-bc-q}
\pd{q}{s}(0,t)= \pd{q}{s}(L,t)=0, \quad t>0.
\end{align}

When the interface is a closed loop, periodic boundary conditions on the node indices and their arc length positions in the discrete model in Eqs (<ref>) translate as periodic boundary conditions on $q(s,t)$:

\begin{align}\label{pbc-q}
q(0,t) = q(L,t), \quad t>0.
\end{align}

Hookean restoring force In Murphy et al., 2019, numerical simulations of the discrete models were performed by rescaling the parameters $a, k, \eta$ in Eq. (<ref>) with $m$. For clarity, we now add a superscript `$(m)$' to these quantities. The scalings were defined as follows:

\begin{align}
&a^{(m)} =\frac{a^\ast}{m}, \qquad k^{(m)} = k^\ast m,
\qquad \eta^{(m)} = \frac{\eta^\ast}{m}, \label{rescale-a-k-eta}
\end{align}

where $a^\ast$, $k^\ast$, and $\eta^\ast$ are $\Order(1)$ as $m\to\infty$ and can be considered mechanical properties of the cell as a whole. Similar parameter rescalings are also mentioned in Murray et al. [Murray et al., 2009] without justification. These rescalings ensure that simulations of the discrete model with linear restoring force match simulations of the continuum cell density. The scaling of the resting spring length $a^{(m)}$ is suggested by springs being smaller with increasing $m$. The scaling of the spring constant $k^{(m)}$ is suggested by $k^\ast=k^{(m)}/m$ being the equivalent spring constant for $m$ identical Hookean springs in series. The reduction of the drag coefficient $\eta^{(m)}$ with $m$ can be argued based on needing to retain a constant total drag force on the cell as its number of springs increases. Equation (<ref>) provides a mathematical justification for these scalings as they ensure that the limit is well defined. Indeed, with the Hookean restoring force law in Eq. (<ref>), we have

\begin{align}
F(q(s,t)) = \lim_{m\to\infty}\frac{k^{(m)}}{m^2\eta^{(m)}}\left(\frac{1}{q(s,t)} - ma^{(m)}\right)= \frac{k^\ast}{\eta^\ast}\left(\frac{1}{q(s,t)} - a^\ast\right).\label{F}
\end{align}

The limit in Eq. (<ref>) defines $\eta^\ast F$ as the restoring force law of a single cell with cell resting length $a^\ast$ and cell spring constant $k^\ast$, where $\eta^\ast$ is the drag coefficient of the cell. Substituting the expression for $F(q)$ obtained from Eq. (<ref>) into Eq. (<ref>) shows that cell density relaxes mechanically along the curved interface according to the nonlinear diffusion equation

\begin{align}\label{evo-cell-density-Hooke}
\pd{q}{t} = \pd{}{s}\left(D(q)\pd{q}{s}\right), \qquad D(q)=\frac{k^\ast}{\eta^\ast}\frac{1}{q^2}.
\end{align}

Equation (<ref>) is the same as that found for cells arranged along the $x$-axis except it involves the arc length coordinate $s$ instead of $x$, see Eq. (<ref>) [Murray et al., 2009, Murray et al., 2012, Murphy et al., 2019]. We also note that with the scalings in Eqs (<ref>), the limit in Eq. (<ref>) is reached for any value of $m$, so the order of approximation of Eq. (<ref>) is indeed $\Order(1/m^2)$.

Tangential stress $\sigma_{\tau\tau}$ in the cell body, given by Eq. (<ref>), also converges to a well-defined limit as $m\to\infty$, since the length $\ell$ of a spring scales as $\Order(m^{-1})$:

\begin{align}
\sigma_{\tau\tau} = E\left(\frac{\ell}{a^{(m)}}-1\right) \longrightarrow E\left(\frac{1/q}{a^\ast}-1\right), \qquad m\to\infty.
\end{align}

Nonlinear restoring forces
With the nonlinear restoring force (<ref>) and the parameter scalings in Eqs (<ref>), the limit (<ref>) becomes

\begin{align}\label{nonlinear-restoring-force-cell}
F\big(q(s,t)\big) = \lim_{m\to\infty}\frac{1}{m\eta^{(m)}} k^{(m)}\big(a^{(m)}\big)^2 \big(\rho_0-mq(s,t)\big) = D_0\big(q_0-q(s,t)\big),
\end{align}

where $D_0=k^\ast (a^\ast)^2/\eta^\ast$ is a diffusion constant, and $q_0=\rho_0/m=1/a^\ast$ is the resting cell density. Substituting the expression for $F(q)$ from Eq. (<ref>) into Eq. (<ref>) shows that cell density relaxes mechanically according to linear diffusion along the curved interface:

\begin{align}
\pd{q}{t} = D_0\pd{^2q}{s^2},
\label{evo-cell-density-lineardiff}
\end{align}

where $\p^2/\p s^2$ corresponds to the one-dimensional Laplace–Beltrami operator along the curve $\b r(s)$ [Berger, 2003]. With the nonlinear restoring force (<ref>), tangential stress in the cell body likewise converges to a finite value as $m\to\infty$ since spring length $\ell$ is $\Order(m^{-1})$:

\begin{align}\label{stress}
\sigma_{\tau\tau} = E\left(1 - \frac{a^{(m)}}{\ell}\right) \longrightarrow E\left(1-\frac{a^\ast}{1/q}\right), \qquad m\to\infty.
\end{align}

If the spring restoring force in Eq. (<ref>) is written with an arbitrary scaling factor $\alpha^{(m)}$ such that $f(\ell)=\alpha^{(m)}\big(1/a^{(m)}-1/\ell\big)$, then $\alpha^{(m)}$ must scale as $\Order(1/m)$ for the continuum limit in Eq. (<ref>) to be well defined. In this case, the diffusion constant is given by $D_0=\alpha^\ast/\eta^\ast$ where $\alpha^\ast=\lim_{m\to\infty}m\alpha^{(m)}$.

With the nonlinear restoring force (<ref>) and the scalings in Eqs (<ref>), the limit (<ref>) becomes

\begin{align}
F\big(q(s,t)\big) = \lim_{m\to\infty}\frac{1}{m\eta^{(m)}} \frac{k^{(m)}\big(a^{(m)}\big)^3}{2} \big(\rho_0^2-m^2q^2(s,t)\big) = \frac{D_0^\ast}{2}\big(q_0^2-q^2(s,t)\big),
\end{align}

where $D_0^\ast = k^\ast(a^\ast)^3/\eta^\ast$. Cell density therefore relaxes mechanically according to a well-known type of porous medium diffusivity:

\begin{align}
\pd{q}{t} = \pd{}{s}\left(D(q)\pd{q}{s}\right), \qquad  D(q) = D_0^\ast\, q.
\label{evo-cell-density-porousdiff}
\end{align}

Equation (<ref>) is the porous medium equation with exponent one for the density dependence of diffusivity, also known as the Boussinesq equation [Boussinesq, 1904, Vázquez, 2006].
Here too, tangential stress converges to a finite value as $m\to\infty$:

\begin{align}
\sigma_{\tau\tau} = \frac{E}{2} \left(1-\left(\frac{a^{(m)}}{\ell}\right)^2\right) \longrightarrow \frac{E}{2}\left( 1- \left(\frac{a^\ast}{1/q}\right)^2\right), \qquad m\to\infty.
\end{align}

If the spring restoring force in Eq. (<ref>) is written with an arbitrary scaling factor $\beta^{(m)}$ such that $f(\ell) = \beta^{(m)}\big(1/(a^{(m)})^2 - 1/\ell^2\big)$, then $\beta^{(m)}$ must scale as $\Order(1/m^2)$ for the continuum limit in Eq. (<ref>) to be well defined. In this case, $D_0^\ast = 2 \beta^\ast/\eta^\ast$ where $\beta^\ast=\lim_{m\to\infty}m^2\beta^{(m)}$.

Table <ref> summarises the relationship between the restoring force laws assumed for the springs, and the corresponding diffusivity of cell density found in the continuum limit. Other nonlinear spring restoring forces can be devised to give a wanted diffusivity function $D(q)$ in the continuum limit. The procedure to find the spring restoring force is to first integrate $F'(q)=-D(q)$, and set the integration constant such that $F(1/a^\ast)=0$. The cell restoring force $\eta^\ast F(q)$ can then be rescaled to a spring restoring force $f(\ell)$ by making the parameters of $\eta^\ast F(q)$ dependent on $m$ based on substituting $a^\ast=m a^{(m)}$ and on utilising the limit (<ref>). Alternatively, linearising $\eta^\ast F(q)$ about $q_0=1/a^\ast$ and matching it to Hooke's restoring force law of the cell in Eq. (<ref>) will provide scalings for the parameters of $\eta^\ast F(q)$ with $m$ via those known for Hooke's law in Eqs (<ref>).

Relationship between spring restoring force law $f(\ell)$ in the discrete model and diffusivity of cell density $D(q)$ in the continuum limit. All cases assume the spring resting length scaling $a^{(m)} = a^\ast/m$ and drag coefficient scaling $\eta^{(m)} = \eta^\ast/m$ from Eqs (<ref>).

$f(\ell)$     Scaling ($m\to\infty$)     $D(q)$

$k^{(m)}(\ell-a^{(m)})$     $k^{(m)}\sim k^\ast m$     $(k^\ast/\eta^\ast)\ q^{-2}$

$\alpha^{(m)}\big((a^{(m)})^{-1}-\ell^{-1}\big)$     $\alpha^{(m)}\sim \alpha^\ast/m$     $(\alpha^\ast/\eta^\ast)$

$\beta^{(m)}\big((a^{(m)})^{-2}-\ell^{-2}\big)$     $\beta^{(m)}\sim \beta^\ast/m^2$     $(2\beta^\ast/\eta^\ast)\ q$

§.§ Numerical simulations

The partial differential equation (<ref>) with boundary conditions (<ref>) or (<ref>) is discretised in space using the method of lines with central finite difference approximations of both partial derivatives involving $s$, and a uniform discretisation of the arc length parameter in $[0,L]$ with $1000$ intervals of length $\deltaup s = L/1000$ [Lynch, 2005]. The resulting system of ordinary differential equations is then solved by the Tsitouras 5/4 Runge–Kutta method using the  method of the Julia package  [Rackauckas & Nie, 2017]. For more detail on the algorithms used to solve the continuum models, the reader is referred to the  computer code available on  [Buenzli et al., 2024].

Comparison of density and stress state between discrete model simulations (magenta) and continuum model simulations (black) for $N=8$ cells around the unit circle with $m=4$ at times $t=0, 0.02, 0.2, 2$ (curved Hookean springs); $R=1$, $k^\ast=1$, $\eta^\ast=1$, $\Deltaup t=0.001$. The initial condition considers that one cell boundary is displaced along the circle by half a resting cell length $a^\ast=2\pi/8$ such that one cell is initially stretched 50% ($\sigma_{\tau\tau}/E=-0.5$), and its neighbouring cell is compressed 50% ($\sigma_{\tau\tau}/E=0.5$), like in Figure <ref>b.

Figure <ref> compares numerical simulations of the continuum partial differential equation (<ref>) with discrete model simulations performed with $N=8$ cells and $m=4$ curved Hookean springs per cell along a circular interface. The evolution of the continuous cell density in Eq. (<ref>) and the evolution of the discrete cell density in Eq. (<ref>) are independent of interface geometry, so the choice of a circular interface is unimportant, except for using periodic boundary conditions. Figure <ref> shows that there is an excellent match between the models for the cell density profiles and for the tangential stress profiles at all times, even for a total number of springs as low as 32. The continuous solution curves cross the stepwise discrete densities and tangential stress near the midpoint between spring boundaries, which is consistent with our continuum limit derivation suggesting $\Order(1/m^2)$ accuracy there as opposed to $\Order(1/m)$ for other points between the spring boundaries. Since the curvature of the interface does not affect the dynamics, these results are similar to those of [Murphy et al., 2019] obtained for springs confined to the $x$ axis, except that periodic boundary conditions are employed here.

Comparison of cell density relaxation with different restoring forces. The discrete model is solved on a circular interface with curved springs for $N=8$ cells and $m=2,4,8$ springs per cells; $R=1$, $k^\ast=1$, $\eta^\ast=1$, $a^\ast=2\pi/8$, $\Deltaup t = 0.001$. The initial condition is the same for all cases and matches that of Figure <ref> ($t=0$) with one elongated cell adjacent to one compressed cell. Density profile along the interface are shown at time $t=0.05$. (a) Continuum cell density profiles obtained by solving Eq. (<ref>) with $D(q)=k^\ast/(\eta^\ast q^2)$, $D(q)=D_0=k^\ast(a^\ast)^2/\eta^\ast\approx 0.62$, and  $D(q)=D_0^\ast q$, where $D_0^\ast=k^\ast(a^\ast)^3/\eta^\ast\approx 0.48$; (b)–(d) Comparison between discrete and continuum model simulations for (b) Hookean springs leading to $D(q)=k^\ast/(\eta^\ast q^2)$; (b) nonlinear springs leading to linear diffusion $D(q)=D_0$; (c) nonlinear springs leading to porous medium diffusion $D(q)=D_0^\ast q$.

In Figure <ref>, we show how the choice of restoring force law affects the dynamics of mechanical relaxation. Figure <ref>a compares the continuous cell densities profiles $q(s,t)$ at time $t=0.05$ obtained by solving Eq. (<ref>) with the diffusivities obtained in the continuum limit for the three choices of restoring force laws in Eqs (<ref>)–(<ref>) (see Table <ref>). The different diffusivities result in slightly different cell density profiles. Figures <ref>b–d show for each choice of restoring force law how quickly increasing the number of springs per cell makes the discrete cell densities converge to the continuum cell density. While the discrete density at the arc length midpoint of a springs is closest to the continuum cell density, it is clear that any other choice of point along the spring will also converge to the continuum density as $m\to\infty$. It is interesting to point out here that the discrete model simulations can be used as an alternative, robust and conservative discretisation method for a broad class of diffusion equations with linear and nonlinear diffusivities.

§.§ Discrete-model justification of parameter rescaling

The scaling with $m$ of the restoring force law in Eq. (<ref>), and the corresponding scalings with $m$ of its parameters in Eqs (<ref>), were shown in Section <ref> to be required for the continuum limit to define a partial differential evolution equation for cell density. In this section, we justify the scalings (<ref>) by analysing the time scales of mechanical relaxation in the evolution equations (<ref>) of the discrete model. For clarity of the presentation, we omit momentarily the superscripts `$(m)$' in $k$, $\eta$, and $a$ in this section.

Hookean springs with periodic boundaries Let $\b s(t)=\big(s_0(t),\ldots,s_{M-1}(t)\big)$ be the vector of arc length positions of the $M$ distinct nodes. With the convention that $s_M(t)=s_0(t)$ due to the periodic boundary condition, and assuming a linear restoring force (Hookean springs), Eqs (<ref>) can be written in matrix form as

\begin{align}\label{discrete-evo-matrix}
\td{}{t}\b s = \frac{k}{\eta}B \b s,  \qquad B=\pmat{-2 & 1 & 0 & \cdots & 0 & 1
\\ 1 & \ddots & \ddots & \ddots &  & 0
\\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots
\\\vdots & \ddots & \ddots & \ddots & \ddots & 0
\\ 0 &  & \ddots & \ddots & \ddots & 1
\\ 1 & 0 & \cdots & 0 & 1 & -2
\end{align}

The matrix $B$ is an $M \times M$ circulant matrix. It is symmetric and therefore has real eigenvalues $\lambda_p$, $p=0,\ldots,M-1$. Decomposing the initial condition in the basis formed by the eigenvectors $\b v^{(p)}$ associated with $\lambda_p$, we can write $\b s(0)=\sum_{p=0}^{M-1}a_p\b v^{(p)}$. The solution to Eq. (<ref>) is then given by

\begin{align}\label{discrete-sol}
\b s(t) = \sum_{p=0}^{M-1} a_p \exp\left(\frac{k}{\eta}\lambda_p t\right)\b v^{(p)}.
\end{align}

To find the characteristic times of the evolution of $\b s(t)$ we now calculate the eigenvalues $\lambda_p$ by solving the eigenvalue problem $B\b v=\lambda \b v$ with $\b v=(v_0,\ldots v_{M-1})$. Given the matrix $B$ in Eq. (<ref>), the components of $\b v$ must satisfy

\begin{align}\label{eigenvalue-problem}
v_{n+1} - (\lambda+2)v_n + v_{n-1} = 0, \qquad n=0,\ldots, M-1,
\end{align}

with the convention that $v_M=v_0$ and $v_{-1}=v_{M-1}$. For $\lambda=0$ the eigenvector $\b v^{(0)}$ solution of Eq. (<ref>) has components $v^{(0)}_n = n$, $n=0,\ldots, M-1$. In Eq. (<ref>), this corresponds to a steady-state contribution $p=0$ in which all the nodes are evenly spaced. To find the other eigenvalues, we take advantage of the periodicity of the system and expand each component $v_n$ of $\b v$ in the Fourier basis $\e^{\i p\,2\pi n/M}$, $p=0,\ldots,M-1$:

\begin{align}\label{fourier}
v_n = \sum_{p=0}^{M-1} \widehat{v}_p\,\e^{\i p\frac{2\pi n}{M}},
\end{align}

so that $v_n = v_{n+M}$ for all $n$. Substituting Eq. (<ref>) into Eq. (<ref>) gives

\begin{align}
\sum_{p=0}^{M-1} \widehat{v}_p\,\e^{\i p\frac{2\pi n}{M}} \left(\e^{\i p\frac{2\pi}{M}} - (\lambda+2) + \e^{-\i p\frac{2\pi}{M}}\right) = 0, \qquad n=0,\ldots,M-1.
\end{align}

For this equality to be satisfied for all $n$ we must have

\begin{align*}
\e^{\i p\frac{2\pi}{M}} - (\lambda + 2) +\e^{-\i p \frac{2\pi}{M}} = 0,\qquad p=0,\ldots,M-1,
\end{align*}

so that

\begin{align}
\lambda_p = -2 + 2 \cos\left(\frac{2\pi p}{M}\right), \qquad p=0,\ldots,M-1.
\end{align}

The largest eigenvalue is $\lambda_0=0$, as found before. All the other eigenvalues are negative, so that the solution (<ref>) always converges to the steady state

\begin{align}
\b{\overline s} = \lim_{t\to\infty}\b s(t) = a_0\b v^{(0)}
\end{align}

in which all nodes are equally spaced. In the long-time limit, the rate of convergence to steady state is dominated by the largest nonzero eigenvalues $\lambda_1=\lambda_{M-1}<0$ (Figure <ref>). In the continuum limit where the number of springs $M\to\infty$, we have

\begin{align*}
\lambda_1 = \lambda_{M-1} =-2+2\cos\left(\frac{2\pi}{M}\right) \sim - \frac{4\pi^2}{M^2}, \qquad M\to \infty,
\end{align*}

so that

\begin{align}
\b s(t) \sim \b{\overline s}  + \big(a_1\b v^\ast + a_{M-1}\b v^{(M-1)}\big) \exp\left(-4\frac{k}{\eta}\frac{\pi^2}{M^2} t\right), \qquad t\to\infty, M\to\infty.
\end{align}

Eigenvalues $\lambda_p$ of the matrix $B$ corresponding to periodic boundary conditions, $p=0,\ldots,M-1$ (plus signs), and eigenvalues $\lambda_p$ of the matrix $\mathcal{B}$ corresponding to fixed boundary conditions, $p=1,\ldots, M-1$ (cross signs), with $M=50$ springs. The inset shows a close-up view of the first eigenvalues. The zero eigenvalue of the matrix $B$ for periodic boundary conditions corresponds to the steady state, see Eq. (<ref>). For fixed boundary conditions, the steady state is not an eigenvector of the matrix $\mathcal{B}$, see Eq. (<ref>).

It is clear from this long-time asymptotic behaviour that the dynamics of mechanical relaxation becomes independent of $M$ in the limit $M\to\infty$ if

\begin{align}\label{asymptotic-params}
\frac{k}{\eta} = \Order\big(M^2\big), \qquad M\to\infty.
\end{align}

This argument can repeated for other pairs of eigenmodes with faster relaxation rates ($\lambda_2=\lambda_{M-2}$, etc). We conclude that the dynamics of the discrete model becomes independent of $M$ in the continuum limit $M\to\infty$, provided that the parameters $k=k^{(m)}$ and $\eta=\eta^{(m)}$ scale with $M$ according to Eq. (<ref>), which is satisfied by the scalings in Eqs (<ref>). The characteristic time scale of relaxation to steady state in this case is given by

\begin{align}\label{relaxation-time-pbc}
T_\text{relax} \sim \frac{1}{4}\frac{\eta^{(m)}M^2}{k^{(m)}\pi^2} = \frac{1}{4}\frac{\eta^\ast N^2}{k^\ast\pi^2}, \qquad m\to\infty.
\end{align}

Hookean springs with fixed boundaries With fixed boundaries, a formal justification of the scalings in Eq. (<ref>) can also be done. The $M$ spring system is such that the arc length position of the first node is fixed at $s_0(t)=s_0=0$, but is now distinct from the arc length position of the last node held fixed at $s_M(t)=s_M=L$.

Let $\b s(t)=\big(s_1(t),\ldots,s_{M-1}(t)\big)$ be the vector of arc length positions of the $M-1$ interior nodes, with the convention that $s_i(t) < s_{i+1}(t)$ for $t \ge 0$ and $i = 0, \ldots, M-1$. With a linear restoring force (Hookean springs), Eqs (<ref>) can be written in matrix form as

\begin{align}\label{discrete-evo-matrix-fixed}
\td{}{t}\b s = \frac{k}{\eta}\mathcal{B} \b s + \dfrac{k}{\eta}\mathcal{A},  \qquad     \mathcal{B}=\pmat{-2 & 1 & 0 & \cdots  & 0
\\ 1 & \ddots & \ddots & \ddots &  \vdots
\\0 & \ddots & \ddots & \ddots &  0
\\ \vdots & \ddots & \ddots & \ddots &  1
\\ 0 & \cdots & 0 & 1 & -2
\mathcal{A} = \pmat[c]{s_0\\
0 \\
\vdots \\
0 \\
} = \pmat[c]{0\\0\\\vdots \\ 0\\ L}.
\end{align}

Here, $\mathcal{B}$ is an $(M-1)\times (M-1)$ Toeplitz tridiagonal matrix and $\mathcal{A}$ is a constant vector of length $M-1$. Writing $\b s(0)=\sum_{p=1}^{M-1}a_p\b v^{(p)}$ in the basis formed by the eigenvectors $\b v^{(p)}$ associated with the (real) $M-1$ eigenvalues $\lambda_p$ of $\mathcal{B}$, the solution of Eq. (<ref>) is

\begin{align}\label{discrete-sol-fixed}
\b s(t) = \b{\overline s} + \sum_{p=1}^{M-2} a_p \exp\left(\frac{k}{\eta}\lambda_p t\right)\b v^{(p)},
\end{align}

where $\b{ \overline s}$ is a constant vector of length $M-1$ whose $i$th element is given by $i L/M$, for $i=1,\ldots,M-1$, which corresponds to the steady state with equally spaced nodes. In this case the eigenvalues of $\mathcal{B}$ can be written as [Kouachi, 2006]

\begin{align}
\lambda_p = -2 + 2\cos\left(\frac{\pi p}{M}\right), \qquad p=1,\ldots,M-1.
\end{align}

Thus, $\lambda_p < 0$ for all $p=1,\ldots M-1$ (Figure <ref>) and all components of the time-dependent term in Eq. (<ref>) decay to zero as $t \to \infty$.  The rate of approach to the long-time limit $\b{\overline s}$ is dominated by the largest eigenvalue

\begin{align*}
\lambda_1 = -2+2\cos\left(\frac{\pi}{M}\right) \sim - \frac{\pi^2}{M^2}, \qquad M\to \infty,
\end{align*}

so that the decay to the steady state is proportional to $\exp\left(- k \pi^2 t /(\eta M^2)\right)$ as $t \to \infty$ and $M \to \infty$.  For the dynamics of mechanical relaxation to become independent of $M$ as $M \to \infty$ we once again must require  $k/\eta = \Order(M^2)$ as in Eq. (<ref>). The characteristic time scale of relaxation with fixed boundaries is

\begin{align}
T_\text{relax}\sim \frac{\eta^{(m)}M^2}{k^{(m)}\pi^2} = \frac{\eta^\ast N^2}{k^\ast\pi^2},
\end{align}

i.e., the dynamics is four times slower than with periodic boundaries (compare with Eq. (<ref>)).

Nonlinear restoring forces
Similar arguments can be made for the long-time relaxation dynamics with nonlinear restoring forces. In this case, Eqs (<ref>) form a nonlinear system of ordinary differential equations that do not have a closed-form solution. To analyse the long-term dynamics of this nonlinear system, we can linearise Eqs (<ref>) about the mechanical equilibrium state $\b{\overline s}=(\overline s_0,\ldots, \overline s_{M-1})$ in which all the nodes are equally spaced: $\overline s_i = iL/M$, where $L$ is the length of the interface. The deviation to steady state $\xi_i = s_i-\overline s_i$ evolves in the linear approximation regime according to

\begin{align}\label{discrete-evo-nonlinear}
\eta\td{\xi_i}{t} = k\big(\xi_{i+1}-2\xi_i+\xi_{i-1}\big),
\end{align}

where $k$ is now defined as $k=f'(L/M)$. Equation (<ref>) can be recast in the same matrix form as Eq. (<ref>) with periodic boundary conditions, or as Eq. (<ref>) with fixed boundary conditions, leading to the time evolutions in Eqs (<ref>) and (<ref>) for $\xi_i$, respectively. It is clear from the expressions of these time evolutions that if $k>0$ then the state with evenly spaced nodes is a linearly stable steady state, and the dynamics of the discrete model near the steady state is independent of $M$ in the continuum limit $M\to\infty$ provided that the parameters $k$ and $\eta$ scale according to Eq. (<ref>). If $k<0$, then the state with evenly spaced nodes is unstable. This situation occurs for example when the restoring force law is chosen to derive from the Lennard–Jones potential, which accounts for aggregation and negative diffusion for a restricted range of density. This restoring force law was found by Murray et al., 2012 to result in evolutions in which small clusters of varying densities emerge since the state with evenly spaced nodes is unstable. In this situation, the linear approximation in Eq. (<ref>) quickly fails and the continuum limit does not represent the discrete model dynamics well [Murray et al., 2012]. In all the restoring forces we consider in this paper, $k>0$. For the nonlinear restoring force $f(\ell) = \alpha^{(m)} (\rho_0-1/\ell)$ with $\rho_0=1/a^{(m)}$, we have

\begin{align}
k^{(m)} = f'\left(\frac{L}{M}\right) = \alpha^{(m)}\frac{M^2}{L^2},
\end{align}

so that with the scalings (<ref>), the parameter $\alpha^{(m)}$ must scale as

\begin{align}
\alpha^{(m)} = k^{(m)}\frac{L^2}{M^2} = \Order\left(\frac{1}{m}\right).
\end{align}

For the nonlinear restoring force $f(\ell)=\beta^{(m)}(\rho_0^2 - 1/\ell^2)$, we have

\begin{align}
k^{(m)} = f'\left(\frac{L}{M}\right) = \beta^{(m)}\frac{M^3}{L^3},
\end{align}

so that the parameter $\beta^{(m)}$ must scale as

\begin{align}
\beta^{(m)} = k^{(m)}\frac{L^3}{M^3} = \Order\left(\frac{1}{m^2}\right).
\end{align}

These scalings are also summarised in Table <ref>.

§ SURFACE TENSION AND NORMAL STRESS

Confluent cells on the curved substrate $\b r(s)$ experience two stress components: the tangential stress $\sigma_{\tau\tau}$ introduced in Eq. (<ref>) due to mechanical interactions with neighbouring cells, and a normal stress $\sigma_{nn}$ due to normal reaction forces exerted by the substrate onto the cell (Figure <ref>). Tangential stress is due to the tangential inner force

\begin{align}
\b F_\tau(s,t) = - \widetilde{f}\left(\rho(s,t)\right)\b \tau(s)
\end{align}

exerting within the cellular layer (Figure <ref>a), i.e.,

\begin{align}\label{tangential-stress2}
\sigma_{\tau\tau} = \frac{F_\tau^\tau}{A} = \frac{\b \tau\vdot\b F_\tau}{A} = - \frac{\widetilde{f}(\rho)}{A} = -\frac{E}{k}\frac{f(\ell)}{a},
\end{align}

as in Eq. (<ref>), where $F_\tau^\tau=\b\tau\vdot\b F_\tau$ is the tangential component of $\b F_\tau$, $A$ is the cross-sectional area over which this force is exerted, and our sign convention is such that tensile stress is negative and compressive stress is positive.

Tangential and normal stresses. (a) Tangential stress $\sigma_{\tau\tau}$ is defined as the tangential component of the inner force $\b F_\tau$ (blue arrow) divided by the cross-sectional surface area $A$; (b) Normal stress $\sigma_{nn}$ is defined as the normal component of the inner force $\b F_n$ (large blue arrow) divided by the contact surface area $w\Deltaup s$, where $w$ is the cell width in the out-of-plane direction. The normal force $\b F_n$ is the net reaction force exerted by substrate on the cell between the arc length positions $s$ and $s+\Deltaup s$ (small blue arrows). This normal force is induced by the tangential forces on curved portions of the interface only (see text for further detail).

The normal reaction forces exerting on the cellular layer are induced by the tangential forces on curved portions of the interface only (Figure <ref>b). They subsist even in the continuum limit when the normal force $\b F^{(\text{n})}$ in Eq. (<ref>) acting on spring boundary nodes in the straight spring model goes to zero. Normal reaction forces acting on the cell body are also present for curved springs which assume no normal force $\b F^{(\text{n})}$ on spring boundaries. The normal stress $\sigma_{nn}$ is generated by the fact that an elongated cell body senses some of the curvature of the interface, such that the tangential forces exerted at the cell boundaries will have different directions and generate a net force in the normal direction (Figure <ref>b). The normal reaction force $\b F_n$ of the substrate is opposite to this net force, and distributed over the contact area between the cell and the substrate. Assuming a cell width $w$ in the out-of-plane direction, the surface area of the cellular layer in contact with the substrate over an arc length $\Deltaup s$ is $w \Deltaup s$, so that

\begin{align}
\sigma_{nn} = \frac{F_n^n}{w\Deltaup s} = \frac{\b n\vdot\b F_n}{w\Deltaup s} = \frac{\b n\vdot\big(\b F_\tau(s+\Deltaup s, t) - \b F_\tau(s,t)\big)}{w\Deltaup s}.
\end{align}

To assign a local, continuous normal stress $\sigma_{nn}(s,t)$ at the arc length position $s$ of the interface, we take $\Deltaup s\to 0$ and also assume the continuum limit $m\to\infty$. In these limits, the net mechanical force per unit length of interface exerting on a portion of the interface of length $\Deltaup s$ is given by

\begin{align}
\lim_{\Deltaup s\to 0}\frac{\b F_\tau(s+\Deltaup s, t)-\b F_\tau(s, t)}{\Deltaup s} &=  \pd{}{s}\b F_\tau(s,t) = - \pd{}{s}\left\{\widetilde{f}\big(\rho(s,t)\big)\right\}\b \tau(s) + \widetilde{f}\big(\rho(s,t)\big)\kappa(s)\b n(s),
\end{align}

where we have used the fact that the curvature of the interface $\kappa(s)$ at the arc length position $s$ is defined such that $\b\tau'(s) = - \kappa(s)\b n(s)$ [Berger, 2003]. Our curvature sign convention is such that $\kappa <0$ where the subtrate is concave, and $\kappa>0$ where the substrate is convex. Therefore,

\begin{align}\label{normal-stress}
\sigma_{nn} = \gamma(s,t)\kappa(s),
\end{align}

\begin{align}\label{surface-tension}
\gamma(s,t) = \frac{\widetilde{f}\big(\rho(s,t)\big)}{w}
\end{align}

is by definition the surface tension of the cellular layer. Equation (<ref>) is similar to the Young–Laplace equation of surface tension of droplets. Because normal stress $\sigma_{nn}$ is induced by tangential forces, $\sigma_{nn}$ is nonzero only if the tangential forces are nonzero and if the interface has some curvature. For example, if the substrate is convex and cells are stretched, $\kappa>0$ and $\widetilde{f}>0$, so that $\sigma_{nn}>0$, which represents compressive normal stress. In contrast, tangential stress $\sigma_{\tau\tau}$ in Eq. (<ref>) is manifestly independent of the geometry of the interface, since the evolution of cell density governed by Eq. (<ref>) is also independent of the shape and curvature of the interface.

Tangential and normal stresses along the open curve $\widetilde{\b r}(u)=\big(u,R\sin(u)\big)$ (solid black curve) using the curved spring model with $N=4, m=4$. Springs are coloured by the tangential stress $\sigma_{\tau\tau}/E$. The contact interface between cells and the substrate is coloured by the normal stress $\sigma_{nn}/(Eh)$. Simulation parameters are as in Figure <ref>, i.e., $R=0.8$, $k=4$, $\eta=0.25$, except resting length $a$ is varied. (a) $a\approx 0.449$ is such that there is no tangential stress in steady state; (b) $a\approx 0.899$ is doubled compared to (a), resulting in compressive tangential stress in steady state; (c) $a\approx 0.225$ is halved compared to (a), resulting in tensile tangential stress in steady state.

Figure <ref> shows snapshots of the discrete model of Figure <ref> at times $t=0$ and $t=4$ for three different values of spring resting length $a$. Both the tangential stress $\sigma_{\tau\tau}/E$ and the normal stress $\sigma_{nn}/(Eh) = \kappa f(\ell)/(ka)$ are shown. We normalise $\sigma_{nn}$ by Young's modulus $E$ and the height $h$ of the cell layer so that $\sigma_{nn}/(Eh)$ only depends on parameters of the discrete spring model. This normalised normal stress is calculated from Eqs (<ref>)–(<ref>) based on the discrete tangential force ${f}(\ell_i)$ exerting on each spring, and the continuous value of curvature along the interface given by

\begin{align}
\kappa(u) = \frac{y''(u)x'(u) - x''(u)y'(u)}{\Big[\big(x'(u)\big)^2 + \big(y'(u)\big)^2\Big]^{3/2}}
\end{align}

for a parametrisation $\widetilde{\b r}(u)=\big(x(u),y(u)\big)$ of the interface [Sethian, 1999].

In Figure <ref>a at time $t=0$, the elongated cell (blue coils) generates compressive (positive) normal stress where the substrate is convex (orange interface region). Similarly, the compressed cell (orange coils) generates compressive normal stress where the substrate is concave. If the spring resting length is increased (Figure <ref>b) or decreased (Figure <ref>c), both the tangential stress and normal stress are modified. At time $t=4$, all cells are relaxed mechanically, so that tangential stress is homogenous along the interface. There is no tangential stress in Figure <ref>a (green coils) since the resting length is chosen such that the springs are at their resting length in mechanical equilibrium. The absence of tangential stress means that there is no normal stress either (green interface). In Figure <ref>b at $t=4$, mechanical equilibrium is reached with compressive tangential stress (orange coils). This compressive tangential stress induces tensile normal stress where the substrate is convex (blue interface region), and compressive normal stress where the substrate is convex (orange interface region). The opposite is found in Figure <ref>c at $t=4$, where there is tensile tangential stress in mechanical equilibrium (blue coils). Tensile tangential stress generates compressive normal stress where the substrate is convex (blue interface region), and tensile normal stress where the substrate is concave (orange–red interface region). In all cases, there is no normal stress wherever curvature is zero (green interface regions).

§ CONCLUSIONS

The development of discrete mathematical models of collective cell mechanics in confluent epithelial layers and their continuum limit have provided useful relationships between mechanical cell properties and tissue-scale cell diffusion properties [Murray et al., 2009, Murray et al., 2011, Murray et al., 2012, Fozard et al., 2010, Murphy et al., 2019, Murphy et al., 2020, Lorenzi, Murray & Ptashnyk, 2020, Murphy et al., 2021, Baker, Parker & Simpson, 2019, Tambyah et al., 2020]. These models are particularly important in mechanobiology, where mechanical stress on a cell may influence its differentiation, phenotype, and behaviour [Opas, 1989, Weinans & Prendergast, 1996, Nelson et al., 2005, Keefer & Desai, 2011, Ladoux & Mège, 2017, Xi et al., 2019, Nelson, 2022]. While the discrete models are easy to formulate and conceptualise, their continuum limit can provide useful analytical insights of the discrete models, such as the emergence of travelling waves [Murphy et al., 2021]. The mechanics and growth of biological tissues is complex due to the fact that they often involve large deformations and differential creation of new material within the tissue [Gamsjäger et al., 2013, Goriely, 2017, Ambrosi et al., 2019]. Linear elasticity may not apply since biological tissues are often subjected to large stresses, including when cells proliferate and when they experience residual stress in equilibrium. Cell-based mathematical models of tissue growth are helpful in deriving growth laws of the continuum mechanics of biological tissues from first principles.

In the present work, we generalise a simple model of the mechanical interactions of cells arranged along a one-dimensional axis, to curved cells arranged as a chain along an arbitrary parametric curve in two-dimensional space, with the specific aim to understand the influence of curvature for the mechanical relaxation of such cellular tissues. We propose a new derivation of the continuum limit of this discrete model based on expansions in terms of small spring lengths, with well-defined cell densities over finite domains. This new derivation allows us to justify how spring restoring force laws must rescale with the number of springs to obtain consistent dynamics in the limit, and to calculate the order of approximation provided by the linear or nonlinear diffusion equation obtained in the limit. We also provide estimates of the timescale of mechanical relaxation of the tissue based on analysing the discrete models, and justify that numerical simulations of discrete models that assume Lennard–Jones elastic potentials in Murray et al., 2012 are linearly unstable.

Our findings show that curvature only influences the dynamics when the mechanical behaviour of the cells is modelled by straight springs, and when the number of these springs is small enough that they may bridge across curved regions of the interface. However, cells are soft matter, partly fluid and partly solid, and their fluid nature means that they can deform to take on the curved shape of their substrate, which is better represented by the curved spring model. The mechanical behaviour of a single cell is in itself a complex combination of many dynamic processes that involve actin polymerisation, adhesion complexes in the cell membrane, and hydrostatic pressure of the cytoplasm [Moeendarbary & Harris, 2014, Ladoux & Mège, 2017]. While we did not model these intracellular mechanical properties in detail, the straight spring model may represent the mechanical behaviour of actin fibres in the cytoskeleton particularly for transmitting tensile forces. The curved spring model may represent the fact that cytosol-filled cell membranes are curved and can fit substrate curvature, while resisting compressive pressure. To consider these two aspects jointly, more realistic cell body shapes that account for cell thickness would need to be considered.

The only quantity that explicitly depends on curvature in the continuum limit is the normal stress of the cells. Tangential forces within the cellular tissue layer generates a surface tension on curved regions of the interface. Biological cells are known to change their behaviour depending on their stress state, so knowing both the tangential and normal components of the stress may be useful for developing mathematical models of cells that account for such behaviours.

Our mathematical model focuses on the mechanical relaxation of a cell monolayer in two-dimensions on a static, curved interface, and as such, could be extended in several directions in future works. The interaction between dynamic evolutions of the interface with the mechanics of cells is an important consideration for tissue growth. Cells on moving boundaries experience curvature-controlled crowding and spreading, which can be countered by mechanical relaxation [Alias & Buenzli, 2017, Alias & Buenzli, 2019]. Normal reaction forces and surface tension may also participate in the evolution of the interface [Ladoux & Mège, 2017, Bidan et al., 2012, Bidan et al., 2012, Bidan et al., 2016, Fratzl et al., 2022], for example in negative pressure wound therapies that apply negative pressure to the wound area [Huang,etal, 2014, Flegg et al., 2020]. Cellular behaviours such as proliferation, death, directed motion along the interface, and competition with other tissues are important to include for modelling more complex, dynamic tissues, including tissue engineered constructs known to be strongly controlled by geometry [Murphy et al., 2019, Buenzli et al., 2020, Hegarty-Cremer et al., 2021, Fratzl et al., 2022, Karakaya et al., 2022].

§.§ Acknowledgments

This research was supported by the Australian Research Council (DP190102545, DP230100025).

§ CONTINUUM LIMIT

To derive the continuum evolution equation (<ref>) from the discrete evolution equation (<ref>) as $m\to\infty$, we introduce a time-dependent parametrisation $\overline{\b r}(u,t)$ of the interface $\b r(s)$ that tracks the evolving midpoint spring positions for constant parameter values $u$, i.e.,

\begin{align}\label{bar-u}
\overline{\b r}(u_i, t) = \b r\big(\overline{s}_i(t)\big),
\end{align}

where $u_i = i\Deltaup u$, $i=0,\ldots, M$ are time-independent, evenly spaced coordinates in a finite parameter space $u\in[0,U)$, and $\overline{s}_i(t)$ are the arc length coordinates of the spring midpoints (Figure <ref>). As $m\to\infty$, $\Deltaup u=U/M\to 0$ and $\overline{s}_{i-1}(t)\sim\overline{s}_i(t)\sim\overline{s}_{i+1}(t)$. This parametrisation of the interface is such that

\begin{align}\label{midpoint-arclength-2}
\overline{s}_i(t) = \int_0^{u_i}\d u\ g(u,t), \qquad g(u,t) = \norm{\pd{\overline{\b r}}{u}(u,t)}.
\end{align}

Spring density at arc length position $\overline{s}_i(t)$, which corresponds to the coordinate $u_i$, is represented by

\begin{align}
\rho\big(\overline{s}_i(t),t\big) = \frac{1}{g(u_i,t)\Deltaup u}.\label{rho-g}
\end{align}

Using Eq. (<ref>) to express $\overline{s}_{i\pm 1}(t)$ in terms of $u_{i\pm1}=u_i\pm\Deltaup u$ and expanding about $u_i$ as $\Deltaup u\to 0$, we obtain (omitting the time dependence of $\overline{s_i}$ to simplify notation),

\begin{align}\label{midpoint-offsets}
\overline{s}_{i\pm1} = \overline{s}_i  \pm g(u_i,t)\Deltaup u + \pd{g}{u}(u_i,t)\frac{\Deltaup u^2}{2} \pm \pd{^2g}{u^2}(u_i,t)\frac{\Deltaup u^3}{6} + \Order\big(\Deltaup u^4\big).
\end{align}

Time-dependent reparametrisation of the interface. The time-dependent parametrisation $\overline{\b r}(u,t)$ of the interface $\b r(s)$ maps constant, evenly spaced coordinates $u_i$ to the time-dependent spring midpoint positions $\b r\big(\overline{s}_i(t)\big)$ (open circles).

With the expressions (<ref>) and $\ell_i^2=1/\rho^2(\overline{s}_i,t) \sim g^2(u_i,t)\Deltaup u^2$, expanding the right hand side of Eq. (<ref>) about $u_i$ as $\Deltaup u\to 0$ gives

\begin{align}\label{evo-continuum-rhs}
\td{}{t}\rho\big(\overline{s}_i,t\big) = - \frac{1}{\eta} \pd{^2}{s^2}\widetilde{f}\big(\rho(\overline{s}_i,t)\big) - \frac{1}{\eta} \frac{1}{g^2(u_i,t)}\pd{g}{u}(u_i,t)\pd{}{s}\widetilde{f}\big(\rho(\overline{s}_i,t)\big) + \Order\left(\Deltaup u^2\right), \qquad m\to\infty.
\end{align}

In the left hand side of Eq. (<ref>), we have, using the chain rule,

\begin{align}
\td{}{t}\rho(\overline{s}_i,t) &= \pd{\rho}{t}(\overline{s}_i,t) + \pd{\rho}{s}(\overline{s}_i,t) \td{\overline{s}_i}{t},
\end{align}

where from Eqs (<ref>), (<ref>)–(<ref>), and (<ref>),

\begin{align}\label{evo-midpoint}
\td{\overline{s}_i}{t} = \frac{1}{2\eta}\left[\widetilde{f}\big(\rho(\overline{s}_{i+1},t)\big) - \widetilde{f}\big(\rho(\overline{s}_{i-1},t)\big)\right] = \frac{1}{\eta} \pd{\widetilde{f}}{s}\big(\rho(\overline{s}_i,t)\big) g(u_i,t)\Deltaup u + \Order\left(\Deltaup u^3\right).
\end{align}

Since from Eqs (<ref>)–(<ref>) $g(u_i,t)\Deltaup u = 1/\rho(\overline{s}_i,t)$ and $(\p\rho/\p s)/\rho = -(\p g/\p s)/g = -(\p g/\p u)/g^2$, the left hand side of Eq. (<ref>) asymptotically becomes:

\begin{align}\label{evo-continuum-lhs}
\td{}{t}\rho(\overline{s}_i,t) &= \pd{\rho}{t}(\overline{s}_i,t) - \frac{1}{\eta} \frac{1}{g^2(u_i,t)} \pd{g}{u}(u_i,t)\pd{\widetilde{f}}{s}\big(\rho(\overline{s}_i,t)\big) + \Order\left(\Deltaup u^3\right)
\end{align}

Equating Eqs (<ref>) and (<ref>) shows that the second term in the right hand sides of both equations cancel, so that

\begin{align}
\pd{\rho}{t}(\overline{s}_i,t) &= -\frac{1}{\eta}\pd{^2}{s^2}\widetilde{f}\big(\rho(\overline{s}_i,t)\big) \left(1+ \Order\left(\Deltaup u^2\right)\right),\label{evo-spring-density-continuous-appx}
\end{align}

which is the same as Eq. (<ref>). This development shows that this partial differential equation is a second-order accurate representation of the discrete model equation (<ref>) as $m\to\infty$, i.e., corrections are $\Order(\Deltaup u^2) = \Order(1/m^2)$.

If $\overline{s}_i(t)$ is not defined as the midpoint between $s_i(t)$ and $s_{i-1}(t)$, corrections of order $\Order(1/m)$ are obtained for the evolution of cell density. Indeed, defining

\begin{align}
\overline{s}_i(t) = \frac{1+\epsilon}{2} s_i(t) + \frac{1-\epsilon}{2} s_{i-1}(t)
\end{align}

for $-1\leq \epsilon \leq 1$, such that the midpoint is obtained when $\epsilon=0$, the expansion in Eq. (<ref>) remains the same. However, the evolution of $\overline{s}_i(t)$ in Eq. (<ref>) now has $\epsilon \Order(\Deltaup u^2)$ corrections. These corrections modify Equation (<ref>) into

\begin{align}\label{evo-continuum-lhs2}
\td{}{t}\rho(\overline{s}_i,t) &= \pd{\rho}{t}(\overline{s}_i,t) - \frac{1}{\eta} \frac{1}{g^2(u_i,t)} \pd{g}{u}(u_i,t)\pd{\widetilde{f}}{s}\big(\rho(\overline{s}_i,t)\big)\notag
\\&+\frac{\epsilon}{2\eta}\left(\frac{1}{\rho^2(\overline{s}_i,t)}\pd{\rho}{s}(\overline{s}_i,t)\pd{^2}{s^2}\widetilde{f}\big(\rho(\overline{s}_i,t)\big) - \frac{1}{\rho^3(\overline{s}_i,t)}\left(\pd{\rho}{s}(\overline{s}_i,t)\right)^2\pd{}{s}\widetilde{f}\big(\rho(\overline{s}_i,t)\big) \right)\notag
\\&+ \Order\left(\Deltaup u^3\right)
\end{align}

Equating Eqs (<ref>) and (<ref>), substituting $\rho = m q$, and using the definition of $F(q)$ in Eq. (<ref>), one obtains

\begin{align}
\pd{q}{t}(s,t) = &-\frac{1}{\eta}\pd{^2}{s^2}F\big(q(s,t)\big) \notag
\\&- \frac{\epsilon}{2\eta}\left(\frac{1}{q^2(s,t)}\pd{}{s}q(s,t)\pd{^2}{s^2}F\big(q(s,t)\big) - \frac{1}{q^3(s,t)}\left(\pd{q}{s}(s,t)\right)^2\pd{}{s}F\big(q(s,t)\big)\right)\frac{1}{m}\notag
\\&+ \Order\left(\frac{1}{m^2}\right).
\end{align}

[Alias & Buenzli, 2017] Alias MA and Buenzli PR (2017) Modeling the effect of curvature on the collective behavior of cells growing new tissue.
Biophys J 112:193–204. <https://doi.org/10.1016/j.bpj.2016.11.3203>.

[Alias & Buenzli, 2018] Alias MA and Buenzli PR (2018) Osteoblasts infill irregular pores under curvature and porosity controls: a hypothesis-testing analysis of cell behaviours. Biomech Model Mechaniobiol 17:1357–1371. <https://doi.org/10.1007/s10237-018-1031-x>.

[Alias & Buenzli, 2019] Alias MA and Buenzli PR (2019) A level-set method for the evolution of cells and tissue during curvature-controlled growth. Int J Numer Meth Biomed Engng 2019:e3279. <https://doi.org/10.1002/cnm.3279>.

[Ambrosi et al., 2019] Ambrosi D, Ben Amar M, Cyron CJ, DeSimone A, Goriely A, Humphrey JD, Kuhl E (2019) Growth and remodeling of living tissues: perspectives, challenges and opportunities. J Roy Soc Interface 16:20190233. <http://dx.doi.org/10.1098/rsif.2019.0233>.

[Baker, Parker & Simpson, 2019] Baker RE, Parker A, Simpson MJ (2019) A free boundary model of epithelial dynamics. J Theor Biol 481:61–74. <https://doi.org/10.1016/j.jtbi.2018.12.025>.

[Basan et al., 2010] Basan M, Risler T, Joanny J-F, Sastre-Garau X, Prost J (2010) Homeostatic competition drives tumor growth and metastasis nucleation. HFSP J 3:265–272. <https://doi.org/10.2976/1.3086732>.

[Berger, 2003] Berger M (2003) A Panoramic View of Riemannian Geometry (Springer). <https://doi.org/10.1007/978-3-642-1>.

[Bidan et al., 2012] Bidan CM, Kommareddy KP, Rumpler M, Kollmannsberger P, Bréchet YJM, Fratzl P, Dunlop JWC (2012) How linear tension converts to curvature: Geometric control of bone tissue growth. PLOS ONE 7:e36336. <https://doi.org/10.1371/journal.pone.0036336>.

[Bidan et al., 2012] Bidan CM, Kommareddy KP, Rumpler M, Kollmannsberger P, Fratzl P, Dunlop JWC (2012) Geometry as a factor for tissue growth: Toward shape optimization of tissue engineering scaffolds, Adv Healthcare Mater 2:186–194. <https://doi.org/10.1002/adhm.201200159>.

[Bidan et al., 2016] Bidan CM, Kollmannsberger P, Gering V, Ehrig S, Joly P, Petersen A, Vogel V, Fratzl P, Dunlop JWC (2016) Gradual conversion of cellular stress patterns into pre-stressed matrix architecture during in vitro tissue growth. J R Soc Interface 13:20160136. <http://dx.doi.org/10.1098/rsif.2016.0136>.

[Boussinesq, 1904] Boussinesq J (1904) Recherches théoriques sur l'écoulement des nappes d'eau infiltrées dans le sol et sur le débit des sources. J Math Pures Appliquées 10:5–78.

[Buenzli et al., 2020] Buenzli P R, Lanaro M, Wong CS, McLaughlin MP, Allenby MC, Woodruff MA, Simpson MJ (2020) Cell proliferation and migration explain pore bridging dynamics in 3D printed scaffolds of different pore size. Acta Biomaterialia, 114:285–295. <https://doi.org/10.1016/j.actbio.2020.07.010>.

[Buenzli et al., 2024] Buenzli PR, Kuba S, Murphy RJ, Simpson MJ (2024) Buenzli2024_Mechanical_cell_interactions_on_curved_interfaces Github repository. <https://github.com/prbuen/Buenzli2024_Mechanical_cell_interactions_on_curved_interfaces>.

[Browning et al., 2021] Browning A, Maclaren OJ, Buenzli PR, Lanaro M, Allenby MC, Woodruff MA, Simpson MJ (2021). Model-based data analysis of tissue growth in thin 3D printed scaffolds, J Theor Biol 528:110852. <https://doi.org/10.1016/j.jtbi.2021.110852>.

[Callens et al., 2020] Callens SJP, Ulyttendaele RJC, Fratila-Apachitei LE, Zadpoor AA (2020) Substrate curvature as a cue to guide spatiotemporal cell and tissue organization, Biomaterials 232:119739, 1–22. <https://doi.org/10.1016/j.biomaterials.2019.119739>.

[Callens et al., 2023] Callens SJP, Fan D, van Hengel IAJ, Minneboo M, Díaz-Payno, Stevens MM, Fratila-Apachitei LE, Zadpoor AA (2023) Emergent collective organization of bone cells in complex curvature fields. Nat Comm 14:855. <https://doi.org/10.1038/s41467-023-36436-w>.

[Chopard & Droz, 1998] Chopard B and Droz M (1998) Cellular automata modeling of physical systems (Cambridge University Press). <https://doi.org/10.1017/CBO9780511549755>.

[Codling, Plank & Benhamou, 2008] Codling EA, Plank MJ, Benhamou S (2008) Random walk models in biology. J Roy Soc Interface 5:813–834. <https://doi.org/10.1098/rsif.2008.0014>.

[Dunlop et al., 2010] Dunlop JWC, Fischer FD, Gamsjäger E, Fratzl P (2010) A theoretical model for tissue growth in confined geometries. J Mech Phys Solids 58:1073–87. <https://doi.org/10.1016/j.jmps.2010.04.008>.

[Flegg et al., 2020] Flegg JA, Menon SN, Byrne HM, McElwain DLS (2020) A current perspective on wound healing and tumour-induced angiogenesis. Bull Math Biol 82:23. <https://doi.org/10.1007/s11538-020-00696-0>.

[Fozard et al., 2010] Fozard JA, Byrne HM, Jensen OE, King JR (2010) Continuum approximations of individual-based models for epithelial monolayers. Math Med Biol 27:39–74. <https://doi.org/10.1093/imammb/dqp015>.

[Fratzl et al., 2022] Fratzl P, Fischer FD, Zickler GA, Dunlop JWC (2022) On shape forming by contractile filaments in the surface of growing tissues. PNAS Nexus 2:1–9. <https://doi.org/10.1093/pnasnexus/pgac292>.

[Gamsjäger et al., 2013] Gamsjäger E, Bidan CM, Fischer FD, Fratyl P, Dunlop JWC (2013) Modelling the role of surface stress on the kinetics of tissue growth in confined geometries. Acta Biomaterialia 9:5531–5543. <http://dx.doi.org/10.1016/j.actbio.2012.10.020>.

[Goriely, 2017] Goriely A (2017) The Mathematics and Mechanics of Biological Growth (Springer). <https://doi.org/10.1007/978-0-387-87710-5>.

[Guyot et al., 2015] Guyot Y, Luyten FP, Schrooten J, Papantoniou I, Geris L (2015) A three-dimensional computational fluid dynamics model of shear stress distribution during neotissue growth in a perfusion bioreactor, Biotech and Bioeng 112:2591–2600. <https://doi.org/10.1002/bit.25672>.

[Hegarty-Cremer et al., 2021] Hegarty-Cremer SGD, Simpson MJ, Andersen TL, Buenzli PR (2021) Modelling cell guidance and curvature control in evolving biological tissues. J Theor Biol 520:110658. <https://doi.org/10.1016/j.jtbi.2021.110658>.

[Huang,etal, 2014] Huang C, Leavitt T, Bayer LR, Orgill DP (2014) Effect of negative pressure wound therapy on wound healing. Curr Prob Surgery 51:301–331. <https://doi.org/10.1067/j.cpsurg.2014.04.001>.

[Jensen & Revell, 2023] Jensen OE and Revell CK (2023) Couple stresses and discrete potentials in the vertex model of cellular monolayers. Biomech Model Mechanobiol 22:1465–1486. <https://doi.org/10.1007/s10237-021-01521-w>.

[Rackauckas & Nie, 2017] Rackauckas C and Nie Q (2017) DifferentialEquations.jl — A performant and feature-rich ecosystem for solving differential equations in Julia. J Open Res Software 5:15. <https://doi.org/10.5334/jors.151>. Software package available from <https://github.com/JuliaDiffEq/DifferentialEquations.jl>.

[Karakaya et al., 2022] Karakaya C, van Asten JGM, Ristori T, Sahlgren CM, Loerakker S (2022) Mechano-regulated cell–cell signaling in the context of cardiovascular tissue engineering. Biomech Model Mechaobiol 21:5–54. <https://doi.org/10.1007/s10237-021-01521-w>.

[Keefer & Desai, 2011] Keefer CL and Desai JP (2011) Mechanical phenotyping of stem cells. Theriogenology 75:1426–1430. <https://doi.org/10.1016/j.theriogenology.2010.11.032>.

[Kouachi, 2006] Kouachi S (2006) Eigenvalues and eigenvectors of tridiagonal matrices. Electronic J Lin Alg 15:115–133. <https://eudml.org/doc/126465>.

[Ladoux & Mège, 2017] Ladoux B and Mège R-M (2017) Mechanobiology of collective cell behaviours. Nat Rev Molec Cell Biol 18:743–757. <https://doi.org/10.1038/nrm.2017.98>.

[LDC2 compiler, 2024] LDC2 compiler of the D programming language, version 1.28.0 (DMD v2.098.0, LLVM 11.1.0) available from <https://github.com/ldc-developers/ldc>, and the standard Phobos runtime library, documented at <https://dlang.org/phobos>.

[Li, Buenzli & Simpson, 2022] Li Y, Buenzli PR, Simpson MJ (2022) Interpreting how nonlinear diffusion affects the fate of bistable populations using a discrete modelling framework. Proc Roy Soc A 478:20220013. <https://doi.org/10.1098/rspa.2022.0582>.

[Li et al., 2022] Li Y, Johnston ST, Buenzli PR, van Heijster P, Simpson MJ (2022) Extinction of bistable populations is affected by the shape of their initial spatial distribution Bull Math Biol 84:21. <https://doi.org/10.1007/s11538-021-00974-5>.

[Lorenzi, Murray & Ptashnyk, 2020] Lorenzi T, Murray PJ, Ptashnyk M (2020) From individual-based mechanical models of multicellular systems to free-boundary problems. Interfaces and Free Boundaries 22:205–244. <https://doi.org/10.4171/ifb/439>.

[Luciano et al., 2021] Luciano M, Xue S-L, De Vos WH, Redondo-Moerata L, Surin M, Lafont F, Hannezo E, Gabriele S (2021) Cell monolayers sense curvature by exploiting active mechanics and nuclear mechanoadaptation. Nat Phys 17:1382–1390. <https://doi.org/10.1038/s41567-021-01374-1>.

[Lynch, 2005] Lynch DR (2005) Numerical partial differential equations for environmnetal scientists and engineers: A first practical course (Springer). <https://doi.org/10.1007/b102052>.

[Martin, Burr & Sharkey, 1998] Martin RB, Burr DB, Sharkey NA (1998) Skeletal Tissue Mechanics (Springer). <https://doi.org/10.1007/978-1-4757-2968-9>.

[Moeendarbary & Harris, 2014] Moeendarbary E, Harris AR (2014) Cell mechanics: principles, practices, and prospects. WIREs Syst Biol Med 6:371–388. <https://doi.org/10.1002/wsbm.1275>.

[Murphy et al., 2019] Murphy RJ, Buenzli PR, Baker RE, Simpson MJ (2019) A one-dimensional individual-based mechanical model of cell movement in heterogeneous tissues and its coarse-grained approximation. Proc Roy Soc A 475:20180838. <https://doi.org/10.1098/rspa.2018.0838>.

[Murphy et al., 2020] Murphy RJ, Buenzli PR, Baker RE, Simpson MJ (2020) Mechanical cell competition in heterogenous epithelial tissues. Bull Math Biol 82:130. <https://doi.org/10.1007/s11538-020-00807-x>.

[Murphy et al., 2021] Murphy RJ, Buenzli PR, Tambuyah TA, Thompson EW, Hugo HJ, Baker RE, Simpson MJ (2021) The role of mechanical interactions in EMT. Phys biol 18:046001. <https://doi.org/10.1088/1478-3975/abf425>.

[Murphy et al., 2021] Murphy RJ, Buenzli PR, Baker RE, Sipmson MJ (2021) Travelling waves in a free boundary mechanobiological model of an epithelial tissue. App Math Lett 111:106636. <https://doi.org/10.1016/j.aml.2020.106636>.

[Murray et al., 2009] Murray PJ, Edwards CM, Tindall MJ, Maini PK (2009) From a discrete to a continuum model of cell dynamics in one dimension. Phys. Rev. E 80:031912. <https://doi.org/10.1103/PhysRevE.80.031912>.

[Murray et al., 2010] Murray PJ, Kang J-W, Mirams GR, Shin S-Y, Byrne HM, Maini PK, Cho K-H (2010) Modelling spatially regulated $\betaup$-catenin dynamics and invasion in intestinal crypts. Biophys J 99:716–725. <https://doi.org/10.1016/j.bpj.2010.05.016>.

[Murray et al., 2011] Murray PJ, Walter A, Fletcher AG, Edwards CM, Tindall MJ, Maini PK (2011) Comparing a discrete and continuum model of the intestinal crypt. Phys Biol 8:026011. <https://doi.org/10.1088/1478-3975/8/2/026011>.

[Murray et al., 2012] Murray PJ, Edwards CM, Tindall MJ, Maini PK (2012) Classifying general nonlinear force law in cell-based models via the continuum limit. Phys Rev E 85:021921. <https://doi.org/10.1103/PhysRevE.85.021921>.

[Nelson et al., 2005] Nelson CM, Jean RP, Tan JL, Liu WF, Sniadecki NJ, Spector AA, Chen CS (2005) Emergent patterns of growth controlled by multicellular form and mechanics. Proc Natl Acad Sci 102:11594–11599. <https://doi.org/10.1073.pnas.0502575102>.

[Nelson, 2022] Nelson CM (2022) Mechanical control of cell differentiation: Insights from the early embryo. Annu Rev Biomed Eng 24:307–322. <https://doi.org/10.1146/annurev-bioeng-060418-052527>.

[Newman, 1980] Newman WI (1980) Some exact solutions to a non-linear diffusion problem in population genetics and comubstion. J Theor Biol 85:325–334.

[Olver, 1974] Olver FWJ (1974) Asymptotics and special functions (Academic Press). <https://doi.org/10.1016/C2013-0-11254-8>.

[Opas, 1989] Opas M (1989) Expression of the differentiated phenotype by epithelial cells in vitro is regulated by both biochemistry and mechanics of the substratum. Develop Biol 131:281–293. <https://doi.org/10.1016/S0012-1606(89)80001-6>.

[King & McCabe, 2003] King JR and McCabe PM (2003) On the Fisher–KPP equation with fast nonlinear diffusion. Proc R Soc Lond A 459:2529–2546. <https://doi.org/10.1098/ rspa.2003.1134>.

[Poujade et al., 2007] Poujade M, Grasland-Mongrain E, Hertzog A, Jouanneau J, Chavrier P, Ladoux B, Buguin A, Silberzan P (2007) Collective migration of an epithelial monolayer in response to a model wound. Proc Natl Acad Sci 104:15988–15993. <https://doi.org/10.1073.pnas.0705062104>.

[Roose et al., 2003] Roose T, Netti PA, Munn LL, Boucher Y, Jain RK (2003) Solid stress generated by spheroid growth estimated using a linear poroelastic model. Microvasc Res 66:204–212. <https://doi.org/10.1016/S0026-2862(03)00057-8>.

[Rumpler et al., 2008] Rumpler M, Woesz A, Dunlop JWC, van Dongen JT, Fratzl P (2008) The effect of geometry on three-dimensional tissue growth, J R Soc Interface 5:1173–1180. <https://doi.org/10.1098/rsif.2008.0064>.

[Sethian, 1999] Sethian JA (1999) Level Set Methods and Fast Marching Methods. 2nd Edition (Cambridge University Press). <www.cambridge.org/9780521642040>.

[Simpson et al., 2010] Simpson MJ, Landman KA, Hughes BD, Fernando AE (2010) A model for mesoscale patterns in motile populations. Physical A 389:1412–1424. <https://doi.org/10.1016/j.physa.2009.12.010>.

[Schamberger et al., 2023] Schamberger B et al. (2023) Curvature in biological systems: Its quantification, emergence, and implications across the scales. Adv Mater 35:2206110. <https://doi.org/10.1002/adma.202206110>.

[Tambyah et al., 2020] Tambyah TA, Murphy RJ, Buenzli PR, Simpson MJ (2020) A free boundary mechanobiological model of epithelial tissues. Proc Roy Soc A 476:20200528. <https://doi.org/10.1098/rspa.2020.0528>.

[Tse et al., 2012] Tse JM, Ghhen G, Tyrrell JA, Wilcox-Adelman SA, Boucher Y, Jain RK, Munn LL (2012) Mechanical compression drives cancer cells toward invasive phenotype. Proc Natl Acad Sci 17:911–916. <https://doi.org/10.1073/pnas.1118910109>.

[VandenHeuvel, Buenzli & Simpson, 2024] VandenHeuvel DJ, Buenzli PR, Simpson MJ (2024) Pushing coarse-grained models beyond the continuum limit using equation learning. Proc Roy Soc A 480:20230619. <https://doi.org/10.1098/rspa.2023.0619>.

[Vázquez, 2006] Vázquez JL (2006) The Porous Medium Equation: Mathematical Theory Oxford Mathematical Monographs (Oxford Academic). <https://doi.org/10.1093/acprof:oso/9780198569039.001.0001>.

[Wang et al., 2018] Wang J, Lo K-Y, McCue S, Simpson MJ (2018) The role of initial geometry in experimental models of wound closing, Chem Eng Science 179:221–226. <https://doi.org/10.1016/j.ces.2018.01.004>.

[Weinans & Prendergast, 1996] Weinans H and Prendergast PJ (1996) Tissue adaptation as a dynamical process far from equilibrium. Bone 19:143–149. <https://doi.org/10.1007/BF00523401>.

[Xi et al., 2019] Xi W, Saw TB, Delacour D, Lim CT, Ladoux B (2019) Material approaches to active tissue mechanics. Nat Rev Mat 4:23–44. <https://doi.org/10.1038/s41578-018-0066-z>.

[Zmurchok, Bhaskar, & Edelstein-Keshet, 2018] Zmurchok C, Bhaskar D, Edelstein-Keshet L (2018) Coupling mechanical tension and GTPase signaling to generate cell and tissue dynamics. Phys Biol 15:046004. <https://doi.org/10.1088/1478-3975/aab1c0>.