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 [Reference Harris, Peter, Bellis, Baum, Kabla and Charras12]. 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 [Reference Meineke, Potten and Loeffler19], 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 [Reference Meineke, Potten and Loeffler19] 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.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.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.