ACTIVE REMODELLING OF TISSUES TO DESCRIBE BIPHASIC RHEOLOGICAL RESPONSES

Abstract

Tissues form from collections of cells that interact together mechanically via cell-to-cell adhesion, mediated by transmembrane cell adhesion molecules. Under a sufficiently large amount of induced stress, these tissues can undergo elastic deformation in the direction of tension, where they then elongate without any topological changes, and experience plastic deformation within the tissue. In this work, we present a novel mathematical model describing the deformation of cells, where tissues are elongated in a controlled manner. In doing so, the cells are able to undergo remodelling through elastic and then plastic deformation, in accordance with experimental observation. Our model describes bistable sizes of a cell that actively deform under stress to elongate the cell. In the absence of remodelling, the model reduces to the standard linear interaction model. In the presence of instant remodelling, we provide a bifurcation analysis to describe the existence of the bistable cell sizes. In the case of general remodelling, we show numerically that cells within a tissue may populate both the initial and elongated cell sizes, following a sufficiently large degree of stress.

1 Introduction

Cell-to-cell adhesion is essential for tissue morphogenesis and homeostasis [21, 23]. Chains of cells within tissues are highly dynamic and can resist the abrasive influences of their environment, including mechanical stresses. The mechanical coupling between lateral surfaces of adjacent cells is mediated by transmembrane cell adhesion molecules (CAMs), which cross the intercellular space and attach to the cell surface or cytoskeleton of neighbouring cells. The process of cell–cell adhesion is dynamic, and the regulation of adhesivity between cells is ongoing due to external mechanical stimuli [21]. A canonical example is the adhesion complex formed between the CAM, cadherin, and the actomyosin cytoskeleton in epithelial tissues [11, 14]. Since the cadherin–actomyosin adhesion complex bears the brunt of mechanical stress between cells in epithelia, it accordingly remodels contact zones between adjacent cells [10, 18]. Under tension, cadherin transmits forces through the cytoskeleton and adapts by either strengthening adhesion complexes to withstand higher forces or dissociating and remodelling adhesion complexes when forces become too large [6, 10], a process which we refer to as “active remodelling.” These dynamics are pronounced under stress relaxation and creep experiments. Here, a monolayered tissue is clamped between two callipers and an external stress is applied, which stretches the tissue, (model depicted in Figure 1). Given a sufficiently large degree of induced stress, the tissue may actively deform elastically, without topological rearrangement, and then undergo plastic deformation [2, 12, 20]. Furthermore, understanding the viscoelastic behaviour of tissues under induced stress has direct consequences for understanding the developmental process of tissue growth, where mechanical forces initial tissue growth, elongation and bending to coordinate morphogenesis [13], and also tissue maintenance [8, 25].

Figure 1 Model schematic, depicting a uni-axial chain of $N=5$ cells undergoing constant elongation between times $t=0$ hrs and $t=\tau $ hrs, after which the tissue is allowed to relax, reaching equilibrium as $t \rightarrow \infty $ hrs. The tissue is clamped at the base (at $x_0 = 0$ ), with the top cell ( $x_N$ ) is free (for $t> \tau $ hrs). The labelling for cells junctions $x_i$ for $0\leq i \leq N$ and cells sizes $l_i$ for $1\leq i \leq N$ are depicted. Initially, the cells all occupy the initial rest lengths (see panel a). However, as the tissue is stretched, some cells remodel, elongating beyond the elongated rest length (see panel b). Upon the release of the applied stretching, the tissue eventually reaches equilibrium with cells occupying a combination of the initial and elongated rest lengths (see panel c).

Mathematical modelling has been used to answer fundamental questions of tissue dynamics and to help better understand a tissue’s behaviours [22]. Traditionally, these models rely on the concept of Hooke’s law and use Newton’s laws of motion, along with overdamped motion, to lead to a force balance, resulting in a simple linear spring interaction between cells with a single fixed point [9, 16, 19]. However, there has been recent work to model bistable springs to describe the mechanical transfer of information [5, 15].

There have been efforts to describe the tissues’ multiphasic response to stress, using a spring-pot model, which behaves as both a spring and a dashpot in different situations [1]. More recently, these spring-pot models have been coupled with traditional models of springs and dashpots, in what is referred to as generalised viscoelastic models [3, 17, 26]. Here, the authors are able to capture the time-dependent response of both single cells and epithelial monolayers under relaxation tests and creep experiments [3]. However, these models are described using fractional calculus which makes their analysis, and subsequent adoption, a challenge.

Another approach to describing the multiscaled response of tissues is to use a power-law response in time. In this type of model, there are multiple time scales describing a cell’s response which allow the biphasic response of the cells to be represented [7]. Alternatively, a more realistic cellular model, the vertex model, can be used to represent cell–cell interactions. Vertex models contain highly detailed cell–cell interactions and intracellular mechanics, making them a suitable model for studying tissue rheology. As such, others have used a modified vertex model with a threshold on tension remodelling and continuous strain relaxation to account for cell remodelling [24].

In this paper, we present a simple individual-based model for a tissue as a uni-axial (1D) chain of cells. In Section 2, we describe the model we propose. In Section 3, we consider three special cases. The first case is no cell remodelling, where the model collapses to the simple linear interaction model. The second case is instant remodelling, where we present a bifurcation analysis to show where the model exhibits mono and bistability. The final case is finite remodelling, which requires numerical methods for analysis. In Section 4, we present these numerical analyses, showcasing how the tissue undergoes remodelling for various remodelling parameters. Finally, in Section 5, we discuss the key results of our model, the limitations of our model and possible future extensions.

2 Model

2.1 Tissue model

In this section, we present our individual-based model for a uni-axial (1D) chain of cells. We consider N cells where the cell junctions are numbered $i=0$ (left boundary) to $i=N$ (right boundary) and are at positions $\{x_i\}_{i=0}^N$ , respectively, such that $x_0(t)$ and $x_N(t)$ are the positions of the boundaries as the tissue grows/stretches and $x_i> x_j$ if $i>j$ . We make the following assumptions:

• • cells are differentiated and so do not undergo proliferation;

• • the timescales considered are sufficiently small, so no topological rearrangement occurs.

These assumptions result in the tissue maintaining a fixed topology, supported by experimental observations [12]. Figure 1 shows the model schematic we will consider throughout this paper.

Each cell i ( $1\leq i \leq N$ ) is assumed to be spring-like and able to remodel its cellular properties, such as its rest length, in response to strain. We denote the size of cell i, $l_i = x_i-x_{i-1}$ for $1\leq i \leq N$ . The value $l_i$ is therefore the length of the cell (between cell junctions i and $i-1$ at positions $x_i$ and $x_{i-1}$ , respectively). The instantaneous behaviour of cell i is Hookean [19], since we assume the tissue to be near a homeostatic state. That is, the tension in the ith cell $T_i$ is given by

(2.1) $$ \begin{align} T_i = \bar{k} ( l_{i} - \Gamma_{i}(t)), \quad 1\leq i \leq N, \end{align} $$

where $\bar {k}$ is the spring constant assumed to be identical for each cell and $\Gamma _i$ is the rest length which we will consider to be dynamic and may be different for each i.

We assume that cell junction movement is over-damped with uniform drag [19] and so any net force $F_i$ placed on the ith cell junction is proportional to the velocity of the junction. That is, since $ {F_i} = T_{i+1} - T_{i}$ (for $1\leq i \leq N-1$ ), (2.1) becomes

(2.2) $$ \begin{align} \frac{d x_i}{dt} = k ( l_{i+1} - l_{i} - \Gamma_{i+1}(t) + \Gamma_{i}(t)), \quad 1\leq i \leq N-1, \end{align} $$

with $l_i = x_i-x_{i-1} \geq 0$ and $\Gamma _i$ the dynamic rest length, for which a model is yet to be defined. Here k is the spring constant scaled by drag.

At this time, we would like to represent the full problem in terms of just $\{l_i\}_{i=1}^N$ , instead of $\{x_i\}_{i=1}^{N-1}$ . Subtracting the ith from the $(i-1)$ th, (2.2) yields

where $\boldsymbol {l} = [l_1, l_2, \ldots , l_{N-1}, l_N]^\dagger $ , and the 1D discrete Laplacian operator for node i is $L_i(\boldsymbol {l}) = l_{i+1}-2l_i+l_{i-1}$ .

Without loss of generality, we can hold the frame of reference with $x_0=0$ . However, the end of the 1D cell $x_N(t)$ may change relative to the origin as the tissue is stretched. For $i=1$ , we may have the clamped (fixed boundary) or free (zero force, $T_{j}=0$ where j is a fictitious external cell) boundary conditions. Using (2.2), we have the equations:

$$ \begin{align*} \frac{d l_1}{dt} &= k ( [{l}_2-\Gamma_2] - (2-\phi_0) [{l}_1-\Gamma_1]), \\[.2cm] \frac{d l_N}{dt} &= k ( [{l}_{N-1}-\Gamma_{N-1}] - (2-\phi_N) [{l}_N-\Gamma_N]) + \phi_N x_N'(t), \end{align*} $$

where in both cases, $\phi =[\phi _0,\phi _N]$ are indicator functions. That is, $\phi _0$ and $\phi _N$ are equal to 1 if the boundaries $x_0$ and $x_N$ , respectively, are clamped, otherwise they take the value zero. Finally, the Nth cell experiences an applied force given by $\phi _N x_N'(t)$ .

Including, therefore, the known boundary conditions $x_0(t)$ and $x_N(t)$ , the model for the chain of cells is determined by solving

(2.3) $$ \begin{align} \frac{d \boldsymbol{l}}{dt} = k\mathcal{L}_\phi ( \boldsymbol{l} - \mathbf{\Gamma} (t)) + \mathbf{b}_\phi \end{align} $$

for some initial $\boldsymbol {l}(0) = \boldsymbol {l}_0$ . Here the discrete 1D Laplacian matrix $\mathcal {L}_\phi $ includes the boundary conditions, given in

$$ \begin{align*} \mathcal{L}_\phi[{1,1}]&=\mathcal{L}_\phi[{N,N}]=-(2-\phi_N), \\ \mathcal{L}_\phi[{1,2}]&=\mathcal{L}_\phi[{N,N-1}]=1, \\ \mathcal{L}_\phi[{i,i}]&=-2 \quad\text{for } 2\leq i\leq N-1, \\ \mathcal{L}_\phi[{i,i-1}]&=\mathcal{L}_\phi[{i,i+1}]=1 \quad\text{for } 2\leq i\leq N-1, \end{align*} $$

and all other elements are zero. The boundaries are included in the matrix $\mathbf {b}_\phi $ ; $\mathbf {b}_\phi [N] = x_N'(t)\phi _N$ , whilst all other elements of $\mathbf {b}$ are zero. Finally, we note that there is nothing in this model that would explicitly force $l_i>0$ for all i. This is because we are assuming that deformations are sufficiently small and the cells remain in the linear regime described by (2.1). At the least, this model is only appropriate if $l_i$ stays strictly positive for each cell.

2.2 Model for dynamic rest length

If a cell is pulled apart by some force to a given length $l_i$ , it is assumed that the cell is capable of remodelling in such a way that, once released, the spring modelling the cell and its external interactions have a new rest length. Of course, since we are assuming that remodelling is a response to strain, as the cells contract towards this new rest length, the strain changes and a new rest length is achieved dynamically through remodelling. We assume therefore that for any given cell size $l_i$ , the cells remodel to a configuration in which the rest length $\gamma _i = \gamma (l_i)$ is desirable. If we assume a tissue is initially at its natural length, $\gamma _0$ , then upon stretching over some nonzero time, the tissue will undergo remodelling. Therefore, upon contraction, the tissue will result in a state greater than its original rest length, $\gamma _0$ . Therefore, we assume that $\gamma '(l_i)>0$ . Furthermore, we assume that for very small deformations, the tissue does not undergo significant changes, $\gamma '(l_i)$ is small and the rest length remains close to its minimum value $\gamma = \gamma _0$ . Stretching the cells of the order of some critical length scale $l_i\sim l_c$ , we assume that the rest length may significantly change up to some maximum rest length $\gamma = \gamma _{{m}}> \gamma _0$ . We model this using the sigmoidal function, giving the equation

(2.4) $$ \begin{align} \gamma_i = \gamma(l_i) = \gamma_0 + (\gamma_{{m}} - \gamma_0 ) \frac{l_i^2}{l_i^2 + l_c^2}. \end{align} $$

This nonlinear response enables us to describe bistable rest lengths. We assume that there is some relaxation time $\lambda ^{-1}$ associated with the remodelling process of each cell and, in response to strain, the dynamic rest length $\Gamma _i(t)$ is determined by the following relaxation equation:

$$ \begin{align*} \frac{d\Gamma_{i}}{dt} = -\lambda [ \Gamma_{i} - \gamma_i] = -\lambda [ \Gamma_{i}(t) - \gamma(l_i(t))], \quad 1\leq i \leq N. \end{align*} $$

In matrix form, this is

(2.5) $$ \begin{align} \frac{d {\mathbf{\Gamma}}}\,{dt} = -\lambda ( \mathbf{\Gamma} - \boldsymbol{\gamma}(\boldsymbol{l})), \end{align} $$

where in this context, $\boldsymbol {\gamma }(\boldsymbol {l}) = [\gamma (l_1), \gamma (l_2), \ldots , \gamma (l_{N-1}), \gamma (l_N)]^\dagger $ .

2.3 Dimensionless model summary

In this manuscript, we are interested in understanding the effects of rest length remodelling on the dynamics of a simple 1D tissue after being stretched. The general model assumes that the parameters can take their full range. The general model consists of (2.3) and (2.5).

Nondimensionalising time with respect to the characteristic timescale associated with spring contraction $k^{-1}$ and length with respect to the separation associated with remodelling $l_c$ , we find the new dimensionless model. To avoid notation clutter, we do not use markings to indicate dimensionless quantities, and instead simply continue with this dimensionless model for the remainder of the manuscript. The following equations give the dimensionless description:

(2.6)

and

(2.7)

with $\boldsymbol {\gamma }(\boldsymbol {l}) = [\gamma (l_1), \ldots , \gamma (l_N)]^\dagger $ as in the nondimensionalising equation (2.4), giving

$$ \begin{align*} \gamma(l) = \gamma_0 + (\gamma_{{m}} - \gamma_0) \frac{l^2}{l^2 + 1} \end{align*} $$

and noting that implies the assumption that the rest lengths are all in equilibrium initially. The vector $\mathbf {b}_\phi =[0,0,\ldots ,0,\phi _N x_N'(t)]^\dagger $ is nonzero, when the right boundary is clamped and stretched or compressed. The remaining parameters of the system include the cell sizes associated with normal contracted $\gamma _0$ and stretched $\gamma _{{m}}$ cells relative to the critical length scale $l_c$ (which is therefore nondimensionalised as 1), $\gamma _0<1<\gamma _m$ . The parameter $\lambda $ describes the rate of remodelling relative to the spring contraction rate.

It is useful for us to study the behaviour of this system in two potential limits. The first of these is a well-studied case where remodelling is slow. This situation is captured in the case where $\lambda \rightarrow 0$ . In this case, we expect that the tissue behaves as an elastic tissue where all of the cells interact by means of linear springs with constant parameters. In the second situation, remodelling is fast $\lambda \rightarrow \infty $ . Indeed, here “fast” and “slow” in a dimensional sense mean that the cells are able to remodel faster or slower than cells are able to contract/stretch, respectively.

We will also find it useful to investigate the behaviour of a single cell. If $N=1$ , then we denote , $\boldsymbol {l}=l_1=x$ and $\boldsymbol {l}_0=\bar {x}$ . The model (2.6)–(2.7) only makes sense in the unclamped regime ( $\phi = [0,0]$ ), because if the cell is clamped at $x_N$ , the solution is trivial ( $x(t) = x_N(t)$ ). In the unclamped regime, the model reduces to

(2.8) $$ \begin{align} \frac{d x}{dt} &= -( x - \Gamma), \quad x(0)=\bar{x}\qquad\qquad\qquad \end{align} $$

(2.9) $$ \begin{align} \frac{d\Gamma}{dt} &= -\lambda ( \Gamma - \gamma(x;\gamma_0,\gamma_{{m}})), \quad \Gamma(0) = \gamma(\bar{x}). \end{align} $$

3 Analysis

3.1 No remodelling ( $\boldsymbol {\lambda = 0}$ )

In the case of no remodelling, the rest length between cells ceases to be dynamic. The solution to (2.7) is . That is, each junction is associated with a static rest length. The tissue then evolves according to

$$ \begin{align*} \frac{d \boldsymbol{l}}{dt} = \mathcal{L}_\phi \boldsymbol{l} + \mathbf{\bar{b}}_\phi, \quad \boldsymbol{l}(0) = \boldsymbol{l}_0, \end{align*} $$

where . This is a relatively uninteresting elastic 1D tissue. The solution $\boldsymbol {l}$ can be found analytically as per the equation

(3.1) $$ \begin{align} \boldsymbol{l}(t) = \exp( \mathcal{L}_\phi t) \bigg[ \int_0^t \exp ( -\mathcal{L}_\phi \tau ) \mathbf{b}(\tau)_\phi \, {d} \tau- \int_0^t \exp ( -\mathcal{L}_\phi \tau) \mathcal{L}_\phi\mathbf{\Gamma}_0 \,{d} \tau + \boldsymbol{l}_0 \bigg]. \end{align} $$

The model here behaves in a very straightforward manner. The tissue relaxes into a steady state if $\lim _{t\rightarrow \infty }\mathbf {b}(t)=0$ or $\phi _N=0$ ; the end cell is not clamped and continuously moving. Since $-\mathcal {L}_\phi $ is symmetric positive definite, the eigenvalues are all positive and therefore as $t\rightarrow \infty $ , the integrals inside the square braces of (3.1) are dominant over the constant $\boldsymbol {l}_0$ . Furthermore, the early dynamic behaviour of $\mathbf {b}_\phi (t)$ adds just a cumulative constant to each element of $\boldsymbol {l}$ . By taking the integrals in (3.1) and recalling the definition of $\mathbf {b}_\phi =[0,0,\ldots ,0,\phi _N x_N'(t)]^\dagger $ , it is possible to show that with clamping at $x_N$ , the long-term behaviour is given by

where $\mathbf {1}$ is a vector of ones and $\lVert \cdot \rVert _1$ is the 1-norm. If there is no clamping, then . Most importantly, this steady state is independent of the initial condition $\boldsymbol {l}_0$ and is therefore globally attracting. It is also the case that in this steady state, the cell sizes are all equal to the initial rest lengths compressed or stretched the same length for each in such a way that the total size of the tissue $\lVert l \rVert _1 = x_N(\infty )$ , as described by the boundary condition. This is because the equilibrium state for a series of linear springs is the state in which the stress is uniform, and in this case, we require the extensions of the cells to be the same since we are assuming they all have the same spring constant.

A single free cell in this regime (2.8)–(2.9) has a length x that simply undergoes exponential decay at a dimensionless rate of 1 to its initial rest length, and is too trivial to warrant further discussion.

3.2 Instantaneous remodelling ( $\boldsymbol {\lambda \rightarrow \infty }$ )

In the case where remodelling is rapid with respect to elastic contraction, we note the pseudo-equilibrium in (2.7), $\mathbf {\Gamma } = \boldsymbol {\gamma }(\boldsymbol {l})$ . Therefore, the system amounts to solving the set of differential equations in

$$ \begin{align*} \frac{d \boldsymbol{l}}{dt} = \mathcal{L}_\phi (\boldsymbol{l} - \boldsymbol{\gamma}(\boldsymbol{l})) + \mathbf{b}_\phi, \quad \boldsymbol{l}(0) = \boldsymbol{l}_0. \end{align*} $$

This is a nontrivial problem due to the nonlinearity $\boldsymbol {\gamma }(\boldsymbol {l})$ , and will require numerical investigation in the next section. We can get some understanding of the behaviour for this system by looking at the behaviour of a single free cell (2.8)–(2.9) as $\lambda \rightarrow \infty $ .

In the case $N=1$ and $\phi = [0,0]$ ,

(3.2) $$ \begin{align} \frac{d x}{dt} = - \bigg[ x - \bigg( \gamma_0 + (\gamma_{{m}} - \gamma_0) \frac{x^2}{x^2 + 1} \bigg) \bigg], \quad x(0)=\bar{x}. \end{align} $$

Such a cell will contract or expand to find equilibrium points. We denote an equilibrium point $x=x_S$ . To find all possible values $x_S$ , we set the time derivative to zero giving

(3.3) $$ \begin{align} 0 = x_S - \bigg( \gamma_0 + (\gamma_{{m}} - \gamma_0) \frac{x_S^2}{x_S^2 + 1} \bigg), \end{align} $$

which rearranges to the following:

$$ \begin{align*} \gamma_0 &= x_S^3 -\gamma_{{m}} x_S^2 + x_S. \end{align*} $$

At this point, we can find the closed-form solution of the cubic solving for $x_S$ ; however, this form is complicated and not particularly insightful. Instead, we define

$$ \begin{align*} f(x_S; \gamma_0, \gamma_{{m}}) &= \gamma_0, \\ g(x_S; \gamma_0, \gamma_{{m}}) &= x_S^3 -\gamma_{{m}} x_S^2 + x_S, \end{align*} $$

and observe when these functions intersect. Noting that g is cubic, we see that it can have at most two stationary points, given by $g'(x_S)=0$ which are at the points

$$ \begin{align*} x_S = \frac{\gamma_{{m}} \pm \sqrt{\gamma_{{m}}^2 - 3}}{3}. \end{align*} $$

We therefore see that there are two different forms for g in terms of stationary points depending on whether $\gamma _m$ is greater than or less than $\sqrt {3}$ , as shown by Figure 2.

Figure 2 Different forms of $g(x_S; \gamma _0, \gamma _{{m}})$ , depending on $ \gamma _{{m}}$ . Here, $g(x_S)$ bifurcates from two stationary points when $\gamma _m>\sqrt {3}$ to no stationary points for $\gamma _m<\sqrt {3}$ .

In the case where $\gamma _m>\sqrt {3}$ , we find it is possible for some $\gamma _0$ values to have multiple fixed points. Graphically, this is shown in Figure 3 by equating $g=f$ and indicating the intersection of these curves on the graph.

Figure 3 Fixed points (black dots) show where f (yellow) and g (purple) intersect (colour available online). These points are solutions to (3.3) and are therefore possible equilibrium sizes for a single cell with rapid remodelling ( $\lambda \rightarrow \infty $ ). If either $\gamma _0$ is outside of a particular window associated with the stationary points of g or if g does not have stationary points, there is only a single fixed point of intersection.

We see that for multiple (three) fixed points $x_S$ , we require that not only should $\gamma _m>\sqrt {3}$ but also $\gamma _0$ be given in the interval

(3.4) $$ \begin{align} \gamma_0 \in \bigg[ \text{max}\bigg\{0, g\bigg(\frac{\gamma_{{m}} + \sqrt{\gamma_{{m}}^2-3}}{3}\bigg)\bigg\}, g\bigg(\frac{\gamma_{{m}} - \sqrt{\gamma_{{m}}^2-3}}{3}\bigg) \bigg] = [ \gamma_0^-(\gamma_m), \gamma_0^+(\gamma_m)]. \end{align} $$

Outside of this region, there is only a single stationary point. We note also by (3.4) that the interval cannot contain negative numbers for $\gamma _0$ .

We show the nature of this bifurcation, see Figure 4. In Figure 4(a), we have plotted the one-parameter bifurcation diagram for $\gamma _0$ , where $\gamma _{m}$ is fixed to $\gamma _{m} = 1.85> \sqrt {3}$ . There are two stable branches of equilibria and an unstable branch of equilibria. We can see the presence of bi-stability for a small range of $\gamma _0$ values (the domain for the unstable branch).

Figure 4 Bifurcation diagram of the single cell with rapid remodelling (3.2). In panel (a), the equilibria branches for $x_S$ as a function of $\gamma _0$ are presented, for a fixed value of $\gamma _{m} = 1.85$ . Solid red curves signify a stable branch of equilibria and the black dashed line signifies the unstable branch for a small region of $\gamma _0$ described by (3.4) (colour available online). In panel (b), the $(\gamma _0,\gamma _m)$ parameter space is split into regions containing only a single stable steady state (these points lie outside of the interval (3.4)) and containing two stable steady states (inside the interval (3.4)). The later region of the parameter space is the region between both blue bifurcation branches which meet at a cusp (beyond which only single stable steady states are expected).

The effect of varying $\gamma _{m}$ on the co-existence of two stable steady states in the system as a function of $\gamma _0$ is given in Figure 4(b). We notice that as $\gamma _0$ is increased, the window of $\gamma _m$ values which permit two stable steady states diminishes to a cusp at $(\gamma _0,\gamma _m) = \sqrt {3}/9 (1,9)$ . We have chosen to plot only the positive quadrant in both bifurcation diagrams. This is because only the region satisfying $0<\gamma _0<\gamma _m$ in the parameter space is allowable, due to physically realistic constraints.

In this case, we expect that if $\gamma _m>\sqrt {3}$ and $\gamma _0\in (\gamma _0^-(\gamma _m),\gamma _0^+(\gamma _m))$ , it will be possible for a cell which starts off with a small rest length on the lower branch to be stretched either by force of the cells around it or by external forces. In this case, the cell can be deformed into the stretched state by crossing the unstable steady state (see the black dashed curves in Figure 4(a)) and at this point, it will tend towards a free rest length on the upper stretched branch. If this happens, to return to the original rest length, compression is necessary (perhaps from neighbouring cells). This proposes an interesting prospect when placing $N>1$ of these cells in series.

3.3 Finite remodelling (finite, nonzero $\boldsymbol {\lambda }$ )

In the case of finite remodelling rates $\lambda $ , if the cells are stretched dynamically beyond the unstable separation distance and released, the cells have time to retract, elastically, towards their small rest length before remodelling has a chance to change the cells substantially enough to push them towards the stretched state at rest. We will see that this can produce some interesting behaviour in the next section, where we re-introduce $\lambda $ finite and $N>1$ to see how a tissue responds to a period of rapid stretching and then release to contract to a new stable steady state.

4 Numerical simulations

Here we will explore numerical simulations of a tissue of N cells described by the model (2.6)–(2.7). All numerical results in this section have been compiled using second-order Runge–Kutta on the model equations. We run numerical experiments on our model undergoing tissue stretching and release. We will run two different types of numerical experiments on this model; stretch and compression.

In our stretch experiments, we will begin with a chain of $N+1$ cells, with each cell initially separated by its minimum rest length $\boldsymbol {l}_0 = \gamma _0\mathbf {1}$ . In this way, the tissue length is $\lVert \boldsymbol {l}_0\rVert _1 = N\gamma _0$ and the cell centres are at $x_i(0) = i\gamma _0$ for each $0\leq i \leq N$ . The tissue is then stretched by clamping the bottom of the tissue $x_0(t) = 0$ and increasing the top of the tissue at a constant speed v for a time $\tau $ , until the size of the tissue has increased by a factor of $\alpha $ . During this time, the cells remodel at a rate $\lambda $ . At time $t=\tau $ , the top of the tissue is released and the tissue is allowed to relax.

4.1 Stretching duration, $\boldsymbol {\tau }$ , and tissue remodelling constant, $\boldsymbol {\lambda }$ , affect tissue elongation dynamics

We first investigate how the tissue behaves under various stretching durations, $\tau $ , with various tissue remodelling constants, $\lambda $ . This will give us an understanding as to how the tissue remodels under various circumstances. For the purpose of this work, we will use the model parameters given in Table 1.

Table 1 Model parameters used to investigate how stretching duration, $\tau $ , and tissue remodelling constant, $\lambda $ , affect tissue elongation dynamics under tissue stretching experiments.

To understand the role stretching has on the remodelling of the tissue, we vary both the stretching duration and the tissue remodelling constant. We will track the tissue size, initial tissue size and also the final tissue size, as well as the (dynamic) cell rest lengths. The final tissue size is calculated by considering the number of cells that are below and above the unstable rest length, and multiplying by the appropriate rest length, as the dynamic remodelling process can require a long time to achieve equilibrium.

We can see immediately that the stretching experiments for the case of slow (no) remodelling ( $\lambda = 0$ ) result in a tissue with the same initial and final tissue length, after the applied stretching event, irrespective of the stretching duration (see Figures 5(a)–5(c)). Consulting the associated cell rest length plots in Figures 6(a)–6(c), we can see that the cell’s rest lengths do not vary, due to the fact that here remodelling does not occur, and so the cells act as simple Hookean springs.

Figure 5 Plot of normalised tissue length ( $l_{\text {Tissue}}$ ) with time (t) for various stretching times $\tau $ and remodelling parameters $\lambda $ . Coloured lines correspond to individual cell’s position, dashed black line is the initial tissue size and solid black line is the steady state tissue size (colour available online).

Figure 6 Plot of normalised cell rest length ( $\Gamma _{\text {Cell}}$ ) with time (t) for various stretching times $\tau $ and remodelling parameters $\lambda $ . Coloured lines correspond to individual cell’s normalised rest length. Lower solid black line is the initial stable cell rest length, the upper solid black line is the elongated stable cell rest length and the dotted line is the unstable cell rest length (colour available online).

In comparison, we observe that as we allow the cells to remodel, we can obtain a larger final tissue. Considering the relatively low strength tissue remodelling constant of $\lambda = 0.002$ , we see that we must stretch the tissue for $\tau = 500$ to achieve a larger than initial tissue (Figure 5(f)), whereas for smaller stretching durations, the tissue returns to its original size, despite its ability to undergo remodelling (see Figures 5(d) and 5(e)). This lack of tissue elongation can be understood by seeing how the cell rest length varies. For small stretching durations, we see that no cells are able to remodel beyond the unstable rest length (Figures 6(d) and 6(e)), and so after the stretching event, these tissues return back to their original size. However, for the larger stretching duration, we can observe that some cells remodel beyond the unstable rest length (Figure 5(f)). As this tissue further relaxes, any cell which has remodelled beyond this unstable rest length will eventually remodel to the elongated stable rest length, resulting in a larger tissue. However, due to the relatively low tissue remodelling constant, this process occurs slowly and therefore the final tissue size in Figure 5(f) is smaller than that shown within the plot.

We next consider a relatively mid strength tissue remodelling constant of $\lambda = 0.02$ . Here, we can see that irrespective of the stretching duration, all of the tissues elongate upon being stretched (see Figures 5(g)–5(i)). However, the extent of the elongation is drastically affected. Looking at the cell rest lengths (Figures 6(g)–6(i)), we see that cells remodel sufficiently quickly to respond to the tissue stretching, meaning that following on from the stretching duration, most cells that will remodel beyond the unstable cell rest length have already done so. We also now observe the stable nature of both the initial and elongated stable rest length in Figure 6(i), as there are some cells which remodel beyond the elongated stable rest length within the stretching duration, but return to the elongated stable rest length after stretching.

Then, we consider a relatively high strength tissue remodelling constant of ${\lambda = 0.2}$ . We observe that the stretching duration affects the final tissue size less so than previously (see Figures 5(j)–5(l)). Looking at the cell rest length in Figures 6(j)–6(l), we see that this is due to the fact that the cells readily remodel and have the ability to pull cells up even after the stretching duration has concluded (most evident in Figure 6(l), where a cell is clearly seen to remodel beyond the unstable rest length after time $t=500$ ). However, there is still some variation in the final tissue size.

Finally, we consider the case of fast (instant) remodelling ( $\lambda \rightarrow \infty $ ). As with the case of a relatively high strength tissue remodelling constant, we see that the stretching duration has less of an effect on the final tissue size, though it still contributes to the tissue elongation (see Figures 5(m)–5(o)). This is due to the fact that the cells instantly remodel to fill the elongated tissue space (see Figures 6(m)–6(o)). Interestingly, in comparison to previous cases where we observed cells remodelling beyond the elongated rest length and pulling neighbouring cells up, if we look closely at the initial stable rest length, as a cell remodels beyond its elongated rest length, it also pushes both of its neighbouring cells. This is most easy to observe in Figure 6(o) where we see elongated cells pushing the neighbour above them (more purple/darker in colour) upwards and cells below them (more yellow/bright in colour) downwards.

4.2 Initial rest length, $\boldsymbol {\gamma _0}$ , elongation parameter, $\boldsymbol {\gamma _m}$ , and tissue stretching factor, $\boldsymbol {\alpha }$ , affect extent of tissue elongation

We now investigate how the tissue behaves under stretching for various initial rest lengths, $\gamma _0$ , elongation parameters, $\gamma _m$ , and tissue stretching factors, $\alpha $ . This will give us an understanding as to how the tissue remodels, providing an understanding to the onset of elongation and also the extent of elongation. For the purpose of this work, we will use the model parameters given in Table 2. As we have observed, upon stretching, the tissue may respond in a heterogeneous manner, with cells at both the initial and elongated state. To obtain a tissue-level understanding, we will measure the tissue elongation factor, which is given as the proportion of cells in the elongated state.

Table 2 Model parameters used to investigate how initial rest length, $\gamma _0$ , elongation parameter, $\gamma _m$ , and tissue stretching factor, $\alpha $ , affect extent of tissue elongation under tissue stretching experiments.

The initial rest length, $\gamma _0$ , and the elongation parameter, $\gamma _m$ , control the distance between the two stable rest lengths. Therefore, we consider the three cases of when the stable rest lengths are: large separation where $(\gamma _0, \gamma _m) = (0.1, 2.0)$ (see Figure 7(a)), a standard separation where $(\gamma _0, \gamma _m) = (0.15, 1.85)$ (see Figure 7(b)) and small separation where $(\gamma _0, \gamma _m) = (0.186, 1.750)$ (see Figure 7(c)). These parameter choices are specified to produce stable rest lengths of either a large, standard or small separation, and to ensure the unstable rest length is at the midpoint between the standard and small separation cases. As before, we will stretch the tissue by a factor of $\alpha $ for a duration $\tau $ , where the tissue remodels with constant $\lambda $ .

Figure 7 (a)–(c) Bifurcation diagrams for different values of $\gamma _m$ , depicting cells with large separation, standard separation and small separation, when specifying a given $\gamma _0$ for each. (d)–(o) Plots of steady state tissue elongation factor with remodelling ( $\lambda $ ) and stretching duration ( $\tau $ ) for various $\gamma _0$ , $\gamma _m$ and $\alpha $ .

We consider first when the separation between the stable rest lengths is large (see Figure 7(a)). These results are shown in Figures 7(d), 7(g), 7(j) and 7(m). Here, irrespective of the tissue stretching factor $\alpha $ , we see that there is always a region within the $\tau \text {--}\lambda $ parameter space where the tissue does not elongate. We note that if we choose a particular $\tau $ and increase $\lambda $ , the tissue transitions from a state of no elongation to elongated. This indicates that the lack of elongation occurs within this region since the tissue does not remodel sufficiently quickly and therefore no cells can access the elongated rest length. Interestingly, if we fix the remodelling constant $\lambda $ and increase the stretching duration, we would expect that the tissue elongation factor would also increase monotonically. However, we can see there are instances where this does not occur, but instead the tissue elongation factor peaks as $\lambda $ increases and then proceeds to decrease again. This indicates that there is a balancing to achieve the maximum tissue elongation, between how quickly the tissue remodels and how quickly the tissue is stretched. Lastly, we would also expect the fastest remodelling tissues to elongate most. However, if we consider $\alpha = 0.75$ (Figure 7(j)), we see that the tissue elongates fully for small $\lambda $ and large $\tau $ . Recalling our previous results, we saw that for small $\lambda $ and large $\tau $ (that is, Figure 6(f)), all of the cells elongated together and so the tissue transitioned in unison to being elongated. However, for the same $\tau $ , but increasing $\lambda $ , (that is, Figure 6(o)), we saw that cells elongated in isolation and then pushed their neighbouring cell down, overshooting the elongated rest length in the process. In this case, since the stable rest lengths are far apart, we see that only small remodelling parameters can elongate fully.

We now consider the case when the stable rest lengths have a standard separation (see Figure 7(b)). Results are shown in Figures 7(e), 7(h), 7(k) and 7(n). As before, we observe the similar style dynamics where varying the stretching duration for some remodelling constants results in the tissue elongation factor being nonmonotonic. We also observe that in this regime where the stable rest lengths have a standard separation, there is a larger region in the $\tau \text {--}\lambda $ parameter space where the tissue elongates fully. However, as before, we observe that full elongation occurs first within the large $\tau $ , small $\lambda $ region.

Lastly, we consider the stable rest lengths to have a small separation (see Figure 7(c)), results shown in Figures 7(f), 7(i), 7(l) and 7(o). Here, we see the onset of full tissue elongation occurs very readily, even for small elongation factors. In this regime, we also observe that for a fixed remodelling constant, the tissue elongation factor increases monotonically with the stretching duration, $\tau $ . We also observe that here, the region in the $\tau \text {--}\lambda $ parameter space where the tissue elongates fully has again expanded further for all tissue stretching factors considered here.

5 Discussion

Cells within tissues are held together by cell-to-cell adhesion, with these mechanical couplings being mediated by transmembrane cell adhesion molecules. However, as this mediation is an active process, the distribution of cell-adhesion molecules throughout the cell reacts to stressors acting upon the cell. This in turn leads to the active remodelling of cells to actively elongate in the direction of tension.

In this work, we have presented a mathematical model to understand how these tissue and cells remodel under tension. The model significantly expands upon current work, where cells are unable to elongate in a controlled manner. Rather, we propose a novel model of the cell with bistable rest lengths. In the limit of slow remodelling (with $\lambda \rightarrow 0$ ), we have shown how this model reduces down to the standard Hookean spring model of a cell, which is valid when the cell is under no (minimal) tension. We have also shown using a bifurcation analysis in the limit of fast remodelling (with $\lambda \rightarrow \infty $ ) how we can maintain bistable rest lengths.

Lastly, using numerical simulations, we have shown how the general remodelling behaves under various tissue stretching experiments. We have shown how we are able to elongate a tissue with the cells actively remodelling, resulting in cells populating both the initial rest length and the elongated rest length. We then showed that irrespective of the initial rest length and elongation parameters used, there exists a persistent region within the $\tau \text {--}\lambda $ parameter space where tissue elongation is prohibited. We also showed that there exists a balance between how readily the cells remodel and how quickly/slowly the tissue is stretched, with how elongated the tissue is at steady state. Slow remodelling and slow stretching generally lead to the most elongated tissue. Finally, we showed that for a particular remodelling constant $\lambda $ , it is possible to obtain nonmonotonic behaviour in the tissue elongation factor with various stretching durations $\tau $ . These observations provide an interesting hypothesis that has the possibility to be experimentally verifiable.

We believe that this model will provide useful insights into the remodelling of mechanical cell-to-cell interactions. Further model extensions include generalising the model into two and three dimensions. An extension to higher dimensions may prevent limitations of our approach, such as topological changes throughout the tissue. Such extensions may be performed in a number of different ways, depending on model specificity. For centre-based models, the extension follows naturally, where our biphasic model could be implemented to describe the changing natural separation between neighbouring cells. For more sophisticated many-node models (such as vertex models), this extension may be implemented in a similar way to our approach; however, using connections to a central node and the cells boundary instead (see [4]). Other potential future investigations could also include understanding how the tissue stretching heterogeneities are sensitive to both tissue size stretching under applied force and the minimal applied force required to induce remodelling, as well as including more sophisticated mechanical force interactions between cells, such as incorporating cell stiffness into a cell’s ability to remodel. Calibrating the model to real-life biological data could provide insights into the particular realistic model dynamics. Moreover, verifying the hypothesis predicted by this work could also prove useful in better understanding the biological dynamics behind cell and tissue remodelling. Finally, while individual-based models of tissues provide a detailed cellular-level understanding into tissue dynamics, extending our model in the continuous limit would provide broader tissue-level insight. However, this would require significant effort to formulate the continuous analogue of our remodelling rest length.