The role of topology and mechanics in uniaxially growing cell networks

TL;DR

This paper explores how the topology and mechanics of cell networks influence tissue growth and stability, revealing that network structure determines whether stress or strain drives growth, contributing to tissue robustness.

Contribution

It introduces a network model linking tissue topology and local growth laws to overall tissue response, clarifying the mechanical basis of tissue stability.

Findings

Series networks are stress-driven for robustness.

Parallel networks are strain-driven for stability.

Network topology influences growth law dominance.

Abstract

In biological systems, the growth of cells, tissues, and organs is influenced by mechanical cues. Locally, cell growth leads to a mechanically heterogeneous environment as cells pull and push their neighbors in a cell network. Despite this local heterogeneity, at the tissue level, the cell network is remarkably robust, as it is not easily perturbed by changes in the mechanical environment or the network connectivity. Through a network model, we relate global tissue structure (i.e. the cell network topology) and local growth mechanisms (growth laws) to the overall tissue response. Within this framework, we investigate the two main mechanical growth laws that have been proposed: stress-driven or strain-driven growth. We show that in order to create a robust and stable tissue environment, networks with predominantly series connections are naturally driven by stress-driven growth, whereas networks with predominantly parallel connections are associated with strain-driven growth.

Figures

Equations

$$\dot{L}=K_{\omega }({\alpha }_{\text{ext}}L-l^{\ast }),$$

$$\dot{L}=K_{\tau }(T-T^{\ast }).$$

$${\mathsf{B}}^T\mathbf{T}={\mathbf{T}}_{\mathrm{ext}}$$

$$T_i=h_i({\alpha }_i),h_i^{\prime }({\alpha }_i)>0,h_i(1)=0.$$

$${\dot{L}}_i={\mathcal{G}}_{\tau }(T_i)$$

$${\dot{L}}_i={\mathcal{G}}_{\omega }(l_i),$$

$$\mathsf{B}^{T}=\left(\begin{array}[]{cccccc}\phantom{-}1&\phantom{-}0&-1&-1&\phantom{-}0&\phantom{-}0\\ \phantom{-}0&\phantom{-}1&\phantom{-}0&\phantom{-}0&-1&\phantom{-}0\\ \phantom{-}0&\phantom{-}0&\phantom{-}0&\phantom{-}1&\phantom{-}0&-1\\ \phantom{-}0&\phantom{-}0&\phantom{-}1&\phantom{-}0&\phantom{-}1&\phantom{-}1\end{array}\right)\,.$$

$$\mathsf{B}^{T}\bm{T}=\begin{pmatrix}T_{1}-T_{3}-T_{4}\\ T_{2}-T_{5}\\ T_{4}-T_{6}\\ T_{3}+T_{5}+T_{6}\end{pmatrix}=\begin{pmatrix}0\\ 0\\ 0\\ 0\end{pmatrix}\quad\begin{array}[]{c}\text{at }w_{1}\\ \text{at }w_{2}\\ \text{at }w_{3}\\ \text{at }w_{4}\,.\end{array}$$

$$l_1=w_1,l_2=w_2,l_3=w_4-w_1,l_4=w_3-w_1,l_5=w_4-w_2,l_6=w_4-w_3.$$

$${\mathbf{\xi }}_1=(-1,1,0,-1,1,-1{)}^T,{\mathbf{\xi }}_2=(0,0,-1,1,0,1{)}^T.$$

$$\underset{\text{loop 1}}{\underbrace{{\mathbf{\xi }}_1\cdot \mathbf{l}=-l_1+l_2-l_4+l_5-l_6=0}},\underset{\text{loop 2}}{\underbrace{{\mathbf{\xi }}_2\cdot \mathbf{l}=-l_3+l_4+l_6=0}}.$$

$$\mathbf{w}\approx {\mathbf{w}}^{\ast }+\epsilon \mathbf{\omega },\mathbf{T}\approx {\mathbf{T}}^{\ast }+\epsilon \mathbf{\tau },\mathbf{L}\approx {\mathbf{L}}^{\ast }+\epsilon \mathbf{\Lambda }.$$

$$\underset{\text{strain-driven growth}}{\underbrace{\left(\begin{array}{ccc}{\mathsf{B}}^T& 0& 0\\ {\mathsf{I}}_m& -{\mathsf{H}}_{\omega }\mathsf{B}& {\mathsf{H}}_{\Lambda }\\ 0& {\mathsf{K}}_{\omega }{\mathsf{H}}_{\omega }^{1/2}\mathsf{B}& 0\end{array}\right)\left(\begin{array}{c}\mathbf{\tau }\\ \mathbf{\omega }\\ \mathbf{\Lambda }\end{array}\right)=\left(\begin{array}{c}{\mathbf{\tau }}_{\text{ext}}\\ 0\\ \dot{\mathbf{\Lambda }}\end{array}\right)}}\text{or}\underset{\text{stress-driven growth}}{\underbrace{\left(\begin{array}{ccc}{\mathsf{B}}^T& 0& 0\\ {\mathsf{I}}_m& -{\mathsf{H}}_{\omega }\mathsf{B}& {\mathsf{H}}_{\Lambda }\\ {\mathsf{K}}_{\tau }{\mathsf{H}}_{\omega }^{-1/2}& 0& 0\end{array}\right)\left(\begin{array}{c}\mathbf{\tau }\\ \mathbf{\omega }\\ \mathbf{\Lambda }\end{array}\right)=\left(\begin{array}{c}{\mathbf{\tau }}_{\text{ext}}\\ 0\\ \dot{\mathbf{\Lambda }}\end{array}\right)}}.$$

$$({\mathsf{H}}_{\omega }{)}_{ii}=\frac{h_i^{\prime }(l_i^{\ast }/L_i^{\ast })}{L_i^{\ast }}>0,({\mathsf{H}}_{\Lambda }{)}_{ii}=\frac{h_i^{\prime }(l_i^{\ast }/L_i^{\ast })l_i^{\ast }}{(L_i^{\ast }{)}^2}>0$$

$$({\mathsf{K}}_{\omega }{)}_{ii}=({\mathsf{H}}_{\omega }{)}_{ii}^{1/2}{\frac{\partial ({\text{G}}_{\omega }{)}_i}{\partial l_i}|}_{\epsilon =0},({\mathsf{K}}_{\tau }{)}_{ii}=({\mathsf{H}}_{\omega }{)}_{ii}^{-1/2}{\frac{\partial ({\text{G}}_{\tau }{)}_i}{\partial T_i}|}_{\epsilon =0},$$

$$\mathbf{\omega }$$

$$\mathbf{\tau }$$

$$\mathsf{C}={\mathsf{H}}_{\omega }^{1/2}\mathsf{B}\in {\mathbb{R}}^{m\times n}$$

$$\mathsf{P}=\mathsf{C}({\mathsf{C}}^T\mathsf{C}{)}^{-1}{\mathsf{C}}^T\in {\mathbb{R}}^{m\times m}.$$

$$\dot{\mathbf{\Lambda }}={\mathsf{J}}_{\omega }\mathbf{\Lambda },{\mathsf{J}}_{\omega }={\mathsf{K}}_{\omega }\mathsf{PD}$$

$$\dot{\mathbf{\Lambda }}={\mathsf{J}}_{\tau }\mathbf{\Lambda },{\mathsf{J}}_{\tau }=-{\mathsf{K}}_{\tau }({\mathsf{I}}_m-\mathsf{P})\mathsf{D}.$$

$$\mathsf{D}={\mathsf{H}}_{\omega }^{-1/2}{\mathsf{H}}_{\Lambda }$$

$$\Pi ({\mathcal{P}}_1{\mathcal{P}}_2)=\Pi ({\mathcal{P}}_2{\mathcal{P}}_1).$$

$$\Pi (\mathcal{K}\mathcal{P}\mathsf{D})=\Pi (\mathsf{D}\mathcal{K}\mathcal{P})=\Pi (\mathsf{D}\mathcal{K}{\mathcal{P}}^2)=\Pi ({\mathsf{P}}^T\mathsf{D}\mathcal{K}\mathcal{P})=\Pi ({\mathcal{P}}^T\mathcal{K}\mathsf{D}\mathcal{P}).$$

$${\sigma }_j({\mathsf{X}}^T\mathsf{YX})={\theta }_j{\sigma }_j(\mathsf{Y}),j=1,\dots ,m$$

$$0\ge {\sigma }_{\text{max}}({\mathcal{K}}_{\le 0}\mathcal{P}\mathsf{D})\ge {\sigma }_{\text{max}}({\mathcal{K}}_{\le 0}\mathsf{D}).$$

$$0\geq\sigma_{\mathrm{max}}(\mathsf{J}_{\tau})\geq-\sigma_{\mathrm{min}}\big{(}(\mathsf{K}_{\tau})_{\geq 0}\mathsf{D}\big{)}$$

$${\sigma }_{\max }(\mathcal{K}\mathcal{P}\mathsf{D})={\sigma }_{\max }({\mathcal{P}}^T\mathcal{K}\mathsf{D}\mathcal{P})=\underset{{\Vert x\Vert }_2=1}{\max }{x}^T{\mathcal{P}}^T\mathcal{K}\mathsf{D}\mathcal{P}x=\underset{\genfrac{}{}{0ex}{}{y\in \text{range}(\mathcal{P})}{{\Vert y\Vert }_2=1}}{\max }{y}^T\mathcal{K}\mathsf{D}y.$$

$$\underset{\genfrac{}{}{0ex}{}{y\in \text{range}(\mathcal{P})}{{\Vert y\Vert }_2=1}}{\max }{y}^T\mathcal{K}\mathsf{D}y\ge \underset{\dim V=n}{\min }\underset{\genfrac{}{}{0ex}{}{y\in V}{{\Vert y\Vert }_2=1}}{\max }{y}^T\mathcal{K}\mathsf{D}y={\sigma }_n(\mathcal{K}\mathsf{D}),$$

$${\mathcal{G}}_{\omega }={\mathsf{K}}_{\omega }(\mathbf{l}-{\mathbf{l}}^{\ast })={\mathsf{K}}_{\omega }\mathsf{B}(\mathbf{w}-{\mathbf{w}}^{\ast }).$$

$${\mathcal{G}}_{\tau }={\mathsf{K}}_{\tau }(\mathbf{T}-{\mathbf{T}}^{\ast }).$$

Full text

I Introduction

Many biological tissues take cues from their mechanical environment to regulate growth and, in turn, generate mechanical stresses on their surrounding hu01 . At the tissue level, arteries respond to wall shear stress taber2001biomechanics and skeletal muscle is primarily driven by mechanical stretch zollner2012stretching , which has been successfully described by the theory of morphoelasticity that combines growth and remodeling with large mechanical deformations goriely17 ; kuhl2014growingMatter ; amatar10 ; goriely2019growth . However, many biological tissues exhibit a cellular structure, and vastly different network topologies are observed in nature, e.g. cells in epithelial monolayers are polygonal with hexagons occurring most frequently in planar slices (Fig. 1A, farhadifar2007influence ), whereas plant root cells form networks of cuboids like bricks in a wall (Fig. 1B, band2012growth ). Bespoke mechanical modelling is required to capture such inherently discrete geometric structures as they grow.

From a mechanical point of view, cell growth causes a local volumetric and mass increase, pushing and pulling neighboring cells and thus leading to a local buildup of mechanical stress. Such geometric and mechanical confinement of a cell in turn adapts and regulates cell growth. The mutual feedback mechanism between cell growth and mechanical and other fields is called a growth law, and the tissue evolution is termed growth dynamics. Mechanical growth laws have been studied both for patterning processes in developing tissues odell1981mechanical ; shraiman2005mechanical and regulatory processes in mature tissues taber2001biomechanics , both with discrete and continuum models. Their forms have been either phenomenologically and micro-structurally inspired goktepe2010athlete ; taber2008theoretical or based on thermodynamical arguments ambrosi2007stress . The main classification is whether growth is stress-driven or strain-driven goktepe2010multiscale ; goriely17 .

When connecting strain-driven or stress-driven growing cells into a network structure, growth dynamics at the network level can become extremely complex shraiman2005mechanical ; taber2008theoretical . Growth dynamics is highly dependent on the form of the growth law vago09 . For a given network topology and growth law, it is unclear if growth dynamics will drive the cell network towards an equilibrium state and provide homeostasis cyron14 ; latorre2018mechanobiological , grow without bounds, or shrink to oblivion erlich18 .

Here, we investigate the role of both the growth law and network topology in homeostasis, in particular how the topology of the network impacts the stability of an equilibrium state if the underlying growth law is stress-driven or strain-driven. While most biological cell networks are two- or three-dimensional, we study here, as a first step, a network deforming along a single axis. The advantage of this approach is that the analytical tractability enables concrete connections to be made between topology, the form of growth law, and stability and thus leads to a broad understanding of the relative role of mechanics and topology in the dynamics of growth, albeit in a simplified geometry. We note that similar questions have been explored in the context of solute transport in networks, e.g. how the vessels of a network adapt to optimize transport ronellenfitsch2016global , or how network topology dynamically adapts to improve transport marbach2016pruning . Our approach is different from the homogenization of micromechanical network models, which are typically limited to (near-) periodic lattices chenchiah2014energy ; murisic2015discrete . Instead, we combine the dynamics of growth with the algebraic structure of an inherently discrete and, possibly, disordered network.

II Growth dynamics of a single cell

To introduce the concept of growth dynamics, consider first a single cell in isolation, which deforms only along one axis. Its deformed length is $l$, and its rest length is $L$. The force $T=h(\alpha)$ in the cell along the growth axis is related to the elastic stretch $\alpha=l/L$ by a monotonic non-linear constitutive law $h$ and assumed to be homogeneous. A simple strain-driven growth law, $\dot{L}=K_{\omega}(l-l^{*})$, dictates that the cell grows when its length is different from the equilibrium length $l^{*}$ with a growth constant $K_{\omega}$. Assuming that a constant external force is prescribed on the cell, $T=F_{\mathrm{ext}}$. Hence we have

[TABLE]

where $\alpha_{\text{ext}}=h^{-1}(F_{\mathrm{ext}})$ and the growth dynamics will converge to the equilibrium state $L^{*}=l^{*}/\alpha_{\text{ext}}$ if and only if the growth constant $K_{\omega}$ is negative. In that case, a cell with length shorter than the equilibrium length, $l<l^{*}$, grows ($\dot{L}>0$) until the equilibrium length is reached. In the case of stress-driven growth, a linear growth dynamic reads

[TABLE]

Similar reasoning shows that under an externally prescribed cell length, the growth constant $K_{\tau}$ must be positive for the growth dynamics to converge to an equilibrium state.

III Growth dynamics of a network of cells

While the stability of growth dynamics for an isolated cell is determined by the sign of the growth constant, the situation is more complicated when multiple cells are organized into networks, which impose additional constraints through length and force balances. The question is then to establish if the entire network will exhibit homeostasis, i.e. whether small perturbations around an equilibrium state will lead the tissue back to its equilibrium. We investigate here the relationship between network topology and growth laws in a highly idealized geometry of cuboid cells connected at their walls and only allowed to deform along one axis. We also note the analogy between a network of cells and electrical circuits; laws such as Kirchhoff’s circuit laws have a natural equivalent in the mechanical network.

We consider a network of $m$ cells which deform along the horizontal $x$-axis as shown in Fig. 2B. The network topology is described by a directed graph in which $m$ edges represent cells and $n$ nodes represent vertical, force-bearing walls. The leftmost wall is not (not counting the leftmost wall) strang1993introduction ; grady2010discrete . Let $\mathsf{B}\in\{-1,0,1\}^{m\times n}$ be the network’s oriented incidence matrix ($B_{ij}=+1$ if the edge $j$ starts at vertex $i$, $-1$ if it arrives at vertex $i$, 0 if it does not connect vertex $i$). Not including the leftmost wall in $\mathsf{B}$ ensures that the network is fixed in space, making the graph Laplacian $\mathsf{B}^{T}\mathsf{B}$ regular. The force balance at each vertical wall (shown in bold in Fig. 2) is given by

[TABLE]

where $\bm{T}=(T_{1},\ldots,T_{m})^{T}$ describes the forces of cells along the $x$-axis and $\bm{T}_{\mathrm{ext}}$ describes the externally prescribed forces at the walls. Eq. (3) is equivalent to Kirchhoff’s current law if the currents are replaced by the forces. The current lengths of cells are $\bm{l}=(l_{1},\ldots,l_{m})^{T}$ and their unstressed lengths are $\bm{L}=(L_{1},\ldots,L_{m})^{T}$. Not counting the leftmost wall, cells are separated by $n$ walls with horizontal coordinates $\bm{w}=(w_{1},\ldots,w_{n})^{T}$, where $w_{n}$ is the rightmost wall and we have $\bm{l}=\mathsf{B}\bm{w}$. We assume a constitutive relationship between the elastic stretch $\alpha_{i}=l_{i}/L_{i}$ of a cell and its force response $h_{i}$:

[TABLE]

Typical laws are neo-Hookean model, $h(\alpha)=\mu(\alpha^{2}-\alpha^{-1})$, or the Hencky model $h(\alpha)=\mu\log\alpha$ mihai2017characterize . Growth dynamics is modeled by the evolution of the rest lengths $L_{i}(t)$. We define the stress-driven growth law

[TABLE]

to be a function of only the force $T_{i}$, and the strain-driven growth law takes the form

[TABLE]

see goriely17 . The governing equations of growth dynamics combine the force balance (3), the constitutive relationships (4) and either the stress-driven or strain-driven growth laws.

Boundary conditions are either fixed total length or imposed force for the entire network. We note that the fixed length condition is the discrete equivalent to a Dirichlet (i.e. prescribed displacement) condition in a continuous system. Likewise prescribed force is equivalent to a Neumann (i.e. prescribed traction) condition. If a force $F_{\mathrm{ext}}$ is prescribed on the outermost walls, then $\bm{T}_{\mathrm{ext}}=(0,\ldots,0,F_{\mathrm{ext}})^{T}$. If instead the total length $L_{\mathrm{ext}}$ is prescribed, the dimension of the incidence matrix is reduced to $m\times(n-1)$, as the position of the rightmost coordinate $w_{n}=L_{\mathrm{ext}}$ is not an unknown. In this case $\bm{T}_{\mathrm{ext}}$ is a zero vector. In the case of linear stability analysis, both boundary conditions can be taken into account by only changing the dimensions of the incidence matrix. Thus, without loss of generality for stability analysis, we focus on the force prescribed case.

IV Length and force balance of an example network

To better illustrate how forces and lengths are balanced in our framework, and to develop a better intuition for how cell graphs and cell diagrams relate to each other, we state the incidence matrix and resulting balances for the example network shown in Figure 2C. In this example, there are $n=4$ vertical walls (not counting the leftmost wall) and $m=6$ cells. The transpose of the incidence matrix is

[TABLE]

The force balance at each vertical wall is given by (3). With no externally prescribed forces ($\bm{T}_{\mathrm{ext}}=\bm{0}$), the force balances at walls 1–4 are

[TABLE]

In our example, lengths $\bm{l}$ are related to wall coordinates by $\bm{l}=\tilde{A}\bm{w}$, which has components

[TABLE]

A basis for the null space of $\mathsf{B}^{T}$ is given by the vectors

[TABLE]

As noted before, the components $+1$ and $-1$ indicate whether the direction goes with or against the arrow, hence $\bm{\xi}_{1}$ corresponds to the loop 1 (red in Figure 3) which comprises segments $l_{1}$, $l_{2}$, $l_{4}$, $l_{5}$ and $l_{6}$. Similarly, $\bm{\xi}_{2}$ corresponds to loop 2 (blue in Figure 3) which comprises segments $l_{3}$, $l_{4}$ and $l_{5}$. The number of independent loops is $\dim N(B^{T})=m-n+1=2$. Note that the outer big loop 3 (green in Figure 3) which comprises segments $l_{1}$, $l_{2}$, $l_{3}$ and $l_{5}$ is represented by the linear combination $\bm{\xi}_{1}+\bm{\xi}_{2}=(-1,1,-1,0,1,0)^{T}$ and is therefore not an independent loop. In summary, the length balance according to Kirchhoff’s voltage law is

[TABLE]

V Mechanical stability of cell assemblies

The goal is to identify the relevant parameters for which the growth laws are stable, and to quantify the probability of their stability for these parameters given a randomly assembled network. We assume that an equilibrium state $\bm{w}^{*}$, $\bm{T}^{*}$, $\bm{L}^{*}$ of the governing equations exists, satisfying $(\mathcal{G}_{\omega})_{i}(l_{i}^{*})=0$ in the strain-driven case and $(\mathcal{G}_{\omega})_{i}(l_{i}^{*})=0$ in the stress-driven case. We keep the form of stress-strain relationships $h_{i}$ and the functional form of the growth law $\mathcal{G}_{i}$ general. With a small parameter $\varepsilon$, we expand

We start with a homeostatic equilibrium state ($\bm{w}^{*}$, $\bm{T}^{*}$, $\bm{L}^{*}$ for which $\dot{L}_{i}=0$ for all $i$) and perform a linear stability analysis by expanding

[TABLE]

Using the series expansion (12) in the governing equations (3), (4) and either the stress-driven growth law $\dot{L}_{i}=\mathcal{G}_{\tau}(T_{i})$ or the strain-driven growth law $\dot{L}_{i}=\mathcal{G}_{\omega}(l_{i})$, leads to the following system at linear order in $\varepsilon$:

[TABLE]

where the square diagonal real matrices with only positive entries $\mathsf{H}_{\omega},\,\mathsf{H}_{\Lambda}\in\mathbb{R}^{m\times m}$ result from linearising the constitutive law

[TABLE]

The square diagonal real matrices $\mathsf{K}_{\omega},\mathsf{K}_{\tau},\,\in\mathbb{R}^{m\times m}$ are the growth constants which characterise how fast a cell grows in relation to the others,

[TABLE]

where $(\mathcal{G}_{\omega})_{i}=(\mathcal{G}_{\omega})_{i}(l_{i})$ is a strain-driven growth law and $(\mathcal{G}_{\tau})_{i}=(\mathcal{G}_{\tau})_{i}(T_{i})$ is a stress-driven growth law. $\mathsf{I}_{m}$ is an $m\times m$ identity matrix. In both cases, $\bm{\tau}$ and $\bm{\omega}$ can be eliminated,

[TABLE]

leading to the linear system (13). Note that $\mathsf{B}^{T}\mathsf{H}_{\omega}\mathsf{B}$ is invertible, as explained in Section III. The growth dynamics of either growth law in the neighborhood of the equilibrium state $\bm{w}^{*}$, $\bm{T}^{*}$, $\bm{L}^{*}$ is linearly stable if the sign of the largest non-zero eigenvalue of the respective Jacobian is negative. A particularly helpful feature of the incidence matrix formulation of this mechanical problem is the emergence of orthogonal projectors. We define

[TABLE]

of full rank as well as

[TABLE]

We note that $\mathsf{P}$ is the orthogonal projector $\mathsf{P}$ onto the range of $\mathsf{C}$, satisfying $\mathsf{P}^{2}=\mathsf{P}^{T}=\mathsf{P}$, whereas $\mathsf{I}_{m}-\mathsf{P}$ is the orthogonal projector onto the null space of $\mathsf{C}^{T}$.

With wall displacements $\bm{\omega}$ and forces $\bm{\tau}$ eliminated according to (16) and (17), the growth dynamics of both strain-driven and stress-driven laws in the neighborhood of the equilibrium state $\bm{\omega}^{*}$, $\bm{T}^{*}$, $\bm{L}^{*}$ can be written as

[TABLE]

where $\mathsf{J}_{\omega}$ and $\mathsf{J}_{\tau}$ are the Jacobians of the strain-driven and stress-driven cases, respectively. Here

[TABLE]

is a positive diagonal matrix encoding the constitutive laws. There are three distinct contributions to the linearized dynamics: the diagonal matrices $\mathsf{K}_{\omega}$ and $\mathsf{K}_{\tau}$ depend on the growth law; the positive diagonal matrix $\mathsf{D}$ is obtained from the constitutive law; and $\mathsf{P}$ encodes the network topology. The special form of these Jacobian matrices implies that their eigenvalues are real, hence: The dynamics of both stress-driven and strain-driven growth laws cannot have oscillations.

Next, we discuss under which conditions a growth law produces stable tissue-level dynamics (i.e. the largest non-zero eigenvalue of the respective Jacobian matrix is negative). We identify three important parameters: the number $k$ of positive entries of either $\mathsf{K}_{\omega}$ or $\mathsf{K}_{\tau}$, the number of cells $m$, and the number of walls $n$ (not counting the leftmost wall). Then, independently of the constitutive law or the topology of the network we have: Strain-driven growth is stable if all growth constants are negative ($k=0$). Stress-driven growth is stable if all growth constants are positive ($k=m$).

Similarly, we can identify cases where the dynamics is unstable: Strain-driven growth is unstable if $k+n>m$. Stress-driven growth is unstable if $k+n<m$. This is true for all constitutive laws and network topologies. These two results thus generalize to an arbitrary network the link between growth constant and stability of an isolated cell.

VI Stability of growth dynamics for arbitrary network topologies, constitutive laws and growth laws

We also have the topological constraint $m\geq n$ (the network cannot have holes) and $0\leq k\leq m$ (there are between 0 and $m$ cells with positive growth constants). Therefore, in the $(m,k,n)$ parameter space the sufficient conditions for stability and instability delineate two triangular pyramids, bounded by the plane $k+n=m$, where stability or instability will occur (indicated in red and blue Fig. 4).

We now proceed to provide proofs for the three propositions written in italics in the main text (in the following, they are stated as Propositions 2, 4 and 6). They relate to the non-oscillatory nature of the growth dynamics, and provide necessary conditions for stable and unstable growth dynamics independently of topology, constitutive laws and growth laws.

We consider an orthogonal projector $\mathcal{P}\in\mathbb{R}^{m\times m}$, which is idempotent and symmetric: $\mathcal{P}^{2}=\mathcal{P}=\mathcal{P}^{T}$. Let the real diagonal matrix $\mathcal{K}\in\text{diag}(\mathbb{R}^{m})$ have $k$ positive entries. Further, $\mathsf{D}\in\text{diag}\left(\mathbb{R}_{>0}^{m}\right)$ is a diagonal matrix of positive entries (this follows from the requirement that constitutive laws are monotonic, see (2)). Note that $\mathsf{D}$ is the same matrix that appears in $(4)_{2}$ and $(5)_{2}$, whereas the form of the orthogonal projector $\mathcal{P}$ and growth constants $\mathcal{K}$ vary depending on the growth law in consideration. Let us also denote by $\Pi(\mathcal{K})$ the spectrum of $\mathcal{K}$ (i.e., the set of all eigenvalues of $\mathcal{K}$).

Firstly, our goal is to prove that the growth dynamics of both growth laws cannot have oscillatory dynamics. To this end, we first prove a lemma about the matrix $\mathcal{K}\mathcal{P}\mathsf{D}$ which mirrors the structure of both growth laws.

Lemma 1

The eigenvalues of the nonsymmetric matrix $\mathcal{J}=\mathcal{K}\mathcal{P}\mathsf{D}$ are real.

Proof. It follows from (Horn2013MatrixAnalysis, , Thm. 1.3.22) that if $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ are two square matrices, then

[TABLE]

Hence, since $\mathcal{K}$, $\mathcal{P}$ and $\mathsf{D}$ are square matrices, and using the fact that $\mathcal{P}=\mathcal{P}^{2}=\mathcal{P}^{T}$, and $\mathsf{D}$ and $\mathcal{K}$ commute, we have that

[TABLE]

This shows that $\mathcal{K}\mathcal{P}\mathsf{D}$ and $\mathcal{P}^{T}\mathcal{K}\mathsf{D}\mathcal{P}$ have the same eigenvalues (note that $\mathcal{K}\mathsf{D}$ is a diagonal matrix). Since $\mathcal{P}^{T}\mathcal{K}\mathsf{D}\mathcal{P}$ is symmetric, the eigenvalues of the nonsymmetric matrix $\mathcal{K}\mathcal{P}\mathsf{D}$ are real. $\square$

Proposition 2

The dynamics of both stress-driven and strain-driven growth laws cannot have oscillations.

Proof. For the strain-driven growth law, the Jacobian is given by (4)2, that is $\mathsf{K}_{\omega}\mathsf{P}\mathsf{D}$, where $\mathsf{K}_{\omega}$ describes the growth constants (15)1 and $\mathsf{P}$ is the orthogonal projector as defined in (19). Thus Lemma 1 applies with $\mathcal{K}=\mathsf{K}_{\omega}$ and $\mathcal{P}=\mathsf{P}$. For the stress-driven case, the Jacobian (5)2 is $-\mathsf{K}_{\tau}(\mathsf{I}_{m}-\mathsf{P})\mathsf{D}$ where $\mathsf{K}_{\tau}$ are the growth constants (15)2. Once again Lemma 1 applies with $\mathcal{K}=-\mathsf{K}_{\tau}$ and $\mathcal{P}=\mathsf{I}_{m}-\mathsf{P}$. It follows for both growth laws that the eigenvalues of the Jacobians are real, and no oscillations are possible. $\square$

In what follows, we order the eigenvalues of an $m\times m$ matrix $\mathcal{J}$ with only real eigenvalues in increasing order, i.e., $\sigma_{\text{min}}(\mathcal{J})=\sigma_{1}(\mathcal{J})\leq\sigma_{2}(\mathcal{J})\leq\ldots\leq\sigma_{m}(\mathcal{J})=\sigma_{\text{max}}(\mathcal{J})$. We will need the following corollary to Ostrowski’s theorem (see (Horn2013MatrixAnalysis, , Cor. 4.5.11)).

Corollary 3

*Let $\mathsf{Y}\in\mathbb{R}^{m\times m}$ be symmetric and let $\mathsf{X}\in\mathbb{R}^{m\times m}$. Then *

[TABLE]

for some nonnegative real scalar $\theta_{j}$ such that $\sigma_{1}(\mathsf{X}^{T}\mathsf{X})\leq\theta_{j}\leq\sigma_{m}(\mathsf{X}^{T}\mathsf{X})$.

Proposition 4

Strain-driven growth is stable if all growth constants are negative ($k=0$). Stress-driven growth is stable if all growth constants are positive ($k=m$).

*Proof. * Applying Corollary 3 to $\mathsf{X}=\mathcal{P}$ and $\mathsf{Y}=\mathcal{K}\mathsf{D}$, and using Lemma 1 leads to $\sigma_{j}(\mathcal{K}\mathcal{P}\mathsf{D})=\theta_{j}\sigma_{j}(\mathcal{K}\mathsf{D})$ with $\sigma_{1}(\mathcal{P}^{T}\mathcal{P})\leq\theta_{j}\leq\sigma_{m}(\mathcal{P}^{T}\mathcal{P})$. From the properties of $\mathcal{P}$, i.e. $\mathcal{P}^{2}=\mathcal{P}^{T}=\mathcal{P}$, it follows that $\mathcal{P}$ has eigenvalues 0 or 1, thus $0\leq\theta_{j}\leq 1$ for all $j=1,\ldots,m$. In the case that $\mathcal{K}$ has only non-positive entries, so does $\mathcal{K}\mathsf{D}$ and we denote $\mathcal{K}$ with only non-positive entries as $\mathcal{K}_{\leq 0}$. Then

[TABLE]

For the strain-driven growth law ($\mathcal{K}=\mathsf{K}_{\omega}$, $\mathcal{P}=\mathsf{P}$), it follows from Eq. (26) that the strain-driven growth with non-positive growth constants $\mathsf{K}_{\leq 0}$ leads to a stable equilibrium irrespective of the network topology encoded in $\mathsf{P}.$

For the stress-driven growth law ($\mathcal{K}=-\mathsf{K}_{\tau}$, $\mathcal{P}=\mathsf{I}_{m}-\mathsf{P}$), we note that $\sigma_{\mathrm{max}}(-\mathsf{K}_{\geq 0}\mathsf{D})=-\sigma_{\mathrm{min}}(\mathsf{K}_{\geq 0}\mathsf{D})$. Eq. (26) then becomes

[TABLE]

and it follows from Eq. (27) that the stress-driven growth with non-negative growth constants $\mathsf{K}_{\geq 0}$ leads to a stable equilibrium irrespective of the network topology encoded in $\mathsf{P}.$ $\square$

The following proposition prepares our claims about unstable growth dynamics, which are stated in the subsequent Proposition 6.

Proposition 5

Let $\mathcal{P}\in\mathbb{R}^{m\times m}$ be an orthogonal projector of rank $n\leq m$, $\mathcal{K}\in\mathbb{R}^{m\times m}$ be diagonal with $k$ positive eigenvalues, and $\mathsf{D}\in\mathbb{R}^{m\times m}$ be diagonal with positive diagonal entries. If $k>m-n$ then the largest eigenvalue of $\mathcal{K}\mathcal{P}\mathsf{D}$ is positive.

Proof. Recall from (24) that $\mathcal{K}\mathcal{P}\mathsf{D}$ and $\mathcal{P}^{T}\mathcal{K}\mathsf{D}\mathcal{P}$ have the same eigenvalues. Then

[TABLE]

Since $\dim(\text{range}(\mathcal{P}))=\text{rank}(\mathcal{P})=n$,

[TABLE]

where the last equality follows from the Courant-Fischer min-max theorem for symmetric matrices (Horn2013MatrixAnalysis, , Thm. 2.4.6). Since $n>m-k$, we have that $n\geq m-k+1$ so that $\sigma_{n}(\mathcal{K}\mathsf{D})\geq\sigma_{m-k+1}(\mathcal{K}\mathsf{D})$. But $\mathcal{K}$, as well as $\mathcal{K}\mathsf{D}$ have $k$ positive eigenvalues so $\sigma_{m-k+1}(\mathcal{K}\mathsf{D})>0$. Hence we conclude that $\sigma_{\max}(\mathcal{K}\mathcal{P}\mathsf{D})\geq\sigma_{n}(\mathcal{K}\mathsf{D})>0$. $\square$

Proposition 6

Strain-driven growth is unstable if $k+n>m$. Stress-driven growth is unstable if $k+n<m$.

Proof. The strain-driven case follows from applying Proposition 5 with $\mathcal{K}=\mathsf{K}_{\omega}$ and $\mathcal{P}=\mathsf{P}$. For the stress-driven case, Proposition 5 can be applied to the Jacobian $-\mathsf{K}_{\tau}(\mathsf{I}_{m}-\mathsf{P})\mathsf{D}$ with $\mathcal{K}=-\mathsf{K}_{\tau}$ and $\mathcal{P}=\mathsf{I}_{m}-\mathsf{P}$. Indeed, $\mathsf{I}_{m}-\mathsf{P}$ is the orthogonal projector onto the orthogonal complement of $\text{range}(\mathsf{C})$, which is the null space of $\mathsf{C}^{T}$. So \text{rank}(\mathsf{I}_{m}-\mathsf{P})=\dim\big{(}\text{range}(\mathsf{I}_{m}-\mathsf{P})\big{)}=\dim\big{(}(\text{null}(\mathsf{C}^{T})\big{)}=m-n. Since $-\mathsf{K}_{\tau}$ has $m-k$ positive eigenvalues, it follows from Proposition 5 that if $k<m-n$ then the largest eigenvalue of $-\mathsf{K}_{\tau}(\mathsf{I}_{m}-\mathsf{P})\mathsf{D}$ is positive. $\square$

VII Some cells try to grow, others shrink: When stability depends on specifics

We now turn our attention to the intermediate region between the stable (red) and unstable (blue) parameter regions in Fig. 4. In this case, stability depends on the specific network topology, constitutive law, and growth law. Here, we choose a Hencky-type constitutive law mihai2017characterize $h_{i}=\mu_{i}\log(\alpha_{i})$. Our goal is to evaluate numerically $\mathsf{J}_{\omega}$, $\mathsf{J}_{\tau}$ to extract their eigenvalues for a given network topology $\mathsf{B}$.

To construct the Jacobians, two challenges must be met: First, a network topology $\mathsf{B}$ must be generated. Secondly, using $\mathsf{B}$, an equilibrium state given by $\bm{w}^{*}$, $\bm{T}^{*}$ and $\bm{L}^{*}$ can be generated. These two steps allow us to compute all the ingredients for the Jacobians $\mathsf{J}_{\omega}$ and $\mathsf{J}_{\tau}$. We now provide some details on the steps.

VII.1 Random generation of reference network geometries

The stress-free reference geometry (Fig. 2A) is created by taking a rectangular unit cell as a starting point, and then dividing it in a sequence of horizontal and vertical cuts. Fig. 5A demonstrates the process. Starting with a single cell (a square, cell 1) we decide randomly whether to make a horizontal or vertical cut. From the resulting two cells (red and blue), we pick a cell at random (in this case cell 2) and randomly decide one of two cut operations (in this case, a horizontal cut). The process is repeated until the final network has been reached. Random numbers are drawn from a uniform distribution using the function RandomInteger in Wolfram Mathematica${}^{\text{\textregistered}}$ 11.2.

Fig. 5B shows a number of randomly generated topologies, sorted by their number of cells and walls. Each topology has an overlay of the graph representation. Moving one unit upward on this grid jumps to a network with one more horizontal cut (i.e. series connection). Moving one unit to the right adds a vertical cut (i.e. parallel connection). The path along which cuts from Fig. 5A have been made is shown with black arrows. The algorithm presented here of making sequential vertical and horizontal cuts places a restriction on the tiling geometries, and does not cover all possible tilings. This algorithm is used for the numerical results in the following section and in Fig.  6. However, we note that the results of Section VI and Fig. 4 are not affected by this restriction of tilings, as they do not depend on the explicit form of the incidence matrix.

VII.2 Numerical linear stability analysis for a non-linear constitutive law

Here we demonstrate how for a specific growth law and constitutive law, a randomly generated numerical Jacobian for a stress-driven or strain-driven growth law is obtained. This random generation of Jacobians underlies Fig. 6.

We focus on a Hencky type constitutive law $h_{i}(\alpha_{i})=\mu_{i}\log\alpha_{i}$. For a strain-driven growth law, we choose

[TABLE]

The stress-driven growth law is

[TABLE]

The first key challenge is to obtain a randomly generated equilibrium state given by $\bm{w}^{*}$, $\bm{T}^{*}$ and $\bm{L}^{*}$ for the non-linear Hencky constitutive law. This requires a numerical solution of (3) and (4). By definition of an equilibrium state $\dot{\bm{L}}^{*}=0$, the growth dynamics trivially satisfies $\bm{\mathcal{G}}_{\omega}(\bm{l}^{*})=\bm{\mathcal{G}}_{\tau}(\bm{T}^{*})=0$. The remaining equations are the force balance and constitutive laws,

[TABLE]

To solve (30), we prescribe a growth profile $\bm{L}^{*}\in\mathbb{R}^{m}$ which we obtain from a uniform random distribution. The reference lengths $\bm{L}^{*}$ correspond to the individual cells’ preferred sizes. We also obtain a reduced incidence matrix from a random graph, obtained through a procedure described in Section VII.2. Thus $\bm{w}^{*}$ and $\bm{T}^{*}$ can be obtained through non-linear root finding. Once they are known, we calculate $\mathsf{H}_{\omega}$ and $\mathsf{H}_{\Lambda}$ according to Eq. (14):

[TABLE]

Then the matrices $\mathsf{C}$ and $\mathsf{D}$ can be obtained as described in the main text. We randomly generate the growth constants $\mathsf{K}_{\omega}$, $\mathsf{K}_{\tau}$ drawing from a uniform distribution, ensuring with an extra constraint that $0\leq k\leq m$ of its of their entries are positive. Then $\mathsf{P}$, $\mathsf{J}_{\omega}$ and $\mathsf{J}_{\tau}$ can be obtained by (19), (20) and (21). In summary, we have described how an equilibrium state compatible with the constraints of a Hencky type constitutive law can be constructed. This allows us to compute all the ingredients needed for the Jacobians $\mathsf{J}_{\omega}$ and $\mathsf{J}_{\tau}$, which tells us the stability of growth dynamics in the neighborhood of a randomly generated equilibrium state $\bm{w}^{*}$, $\bm{T}^{*}$, $\bm{L}^{*}$.

VIII How likely is it that an equilibrium state in a randomly generated network leads to stable dynamics?

To quantify how likely it is that an equilibrium state for a given growth law is stable, we use Jacobian matrices constructed from randomly generated parameters as described above. We define the probability of a stable equilibrium as the number of randomly generated Jacobian matrices of which the largest non-zero eigenvalue is negative, divided by the total number of randomly generated Jacobian matrices. Fig. 6A and B show the probabilities of stable strain-driven and stress-driven growth dynamics, respectively, for up to 50 cells. For reference, the unstable regions from Fig. 4 are superimposed. The color of each pixel corresponds to the probability of stable growth calculated from 100 randomly generated matrices. Fig. 6A* and B* show the fixed $m=50$ case in detail. The largest regions of high probability for stable strain-driven growth occur for small values of $n$ and $k$ compared to $m$ (bottom left corner of A*), implying that parallel connections are more likely to produce a stable configuration. Conversely, for stress-driven growth (top right corner of B*), series connections are more likely to produce a stable configuration. Fig. 7 confirms this via a one-dimensional slice of the data of Fig. 6, showing how probabilities of stable growth dynamics change as the number of horizontal cuts, i.e. parallel connections, increases.

While we are not tracking possible topological modifications such as cell division or neighbor exchanges, the results of Fig. 6 allows us to compare the probability of stability before and after cell division. Assuming a growth dynamics model, we can make inferences about whether cell division planes are likely to be parallel or perpendicular to the direction of growth. Indeed, Fig. 6 implies that in order for cell division to remain stable, planes parallel to the direction of growth are associated with strain-driven growth, whereas cell division planes perpendicular to the direction of growth are associated with stable stress-driven growth. The relationship between cell division planes and mechanical stability is summarised in Fig. 8.

IX Discussion

Our analysis suggests the following mechanism in uniaxial tissue growth dynamics: stress-driven growth laws lead to series connections at the level of the whole network while strain-driven growth leads to parallel connections. Our approach relates the topological structure of the system to its growth dynamics.

We have developed our model under the geometric restriction that cells grow and deform in only one dimension. While this is a clear limitation of the model, the advantage of this simplification is an analytical description that allows us to make concrete statements, as for instance presented in Figs. 4 and 6, that connect network topology and the local rules of growth dynamics at the cell level to the global dynamic response. Moreover, there exist biological systems that are reasonably well-approximated by such a system, for instance the structure of skeletal muscle. In muscles, sarcomeres have been modeled as one-dimensional active elements connected in a complex layered 1D network structure that forms myofibrils, muscle fibres, and further layers that build up the muscle tissue caruel2018physics . Furthermore, in heart tissues, the addition of sarcomeres in series and parallel has been associated with different loading conditions, see goktepe2010athlete , suggesting that there are microscopic mechanisms that locally change the network topology in response to different end conditions.

Another interesting application of these ideas can be found in the growth dynamics of plant stems and roots. For instance, the Arabidopsis thaliana plant root is mostly a uni-directionally growing network where cells are primarily connected in series with respect to the growth direction. The growth dynamics is well captured by a standard physiological model by Lockhart, in which uni-directional plant cell growth is described as proportional to a difference between forces due to the cell’s turgor pressure and a biologically encoded yield pressure, with a positive proportionality constant jensen2015multiscale . Lockhart’s law is a stress-driven growth law with a positive growth constant. This model description, when paired with the observation that the Arabidopsis thaliana is mostly comprised of series connections, is consistent with the criteria of stable stress-driven growth dynamics, and thus supports the hypothesis that tissues maximize the probability of stable growth dynamics.

Aside from dimensionality, our model has several limitations that need to be addressed in a more comprehensive theory of the dynamics of growing biological tissue. Friction and/or adhesion between cells sharing a horizontal wall would be a useful addition, as cortical tension and cell adhesion have been recognised as important factors in epithelial monolayers farhadifar2007influence . Also, cell wall anisotropy jensen2015multiscale and intracellular fluxes cheddadi2019coupling are of key importance in plants. A popular class of models for both types of tissue are vertex models (goriely17, , p. 56).

On the other hand, incorporating some of the components described above would add significant complexity to our model system, whereas the relative simplicity and analytical tractability may be helpful to study the qualitative impact of growth and mechanics in complex biological problems, such as size regulation and growth termination, i.e. how the overall domain size of an organ or tissue (e.g. the wing disc of Drosophila melanogaster) is communicated to the individual cells to stop growth and division. Various cell-based shraiman2005mechanical ; hufnagel2007mechanism and continuum aegerter2007model ; ambrosi2015active models, incorporating a combination of mechanical feedback and morphogen diffusion, postulate different local growth laws and illustrate the complexity of finding a minimal mechanism that incorporates size regulation and growth termination. We believe that our simplified one-dimensional model, to which discrete modeling of advection and diffusion can be coupled, is a first step for theoretical studies of size regulation, as finite equilibrium states of growth dynamics can be interpreted as states of growth termination.

In the light of recent work work by Jensen et al. jensen2019force , an extension of our 1D model to 2D vertex model geometries (exploiting their algebraic structure with incidence matrices and discrete calculus) may be promising. Further extensions of our model to better approximate biological tissue (for instance for the study of growth termination in epithelia) would include advection, diffusion and dilution of chemical patterns aguilar2018critical ; erlich2018physicalDeterminants ; grady2010discrete and would take account of the (near-) incompressibility of cells harris2012characterizing , as well as dynamical topological changes caused by cell division xu2015changes , which requires an accompanying topological growth law.

Ethics statement. This work did not involve any active collection of human data.

Data accessibility statement. All data needed to evaluate the results and conclusions are present in the paper. The associated datasets and codes can be accessed via the Figshare repository (https://doi.org/10.6084/m9.figshare.10315775).

Competing interests statement. We have no competing interests.

Authors’ contributions. AG, DEM, and AE devised the study, AE conducted the analysis, GWJ and AE devised the network description, FT and AE formulated the algebraic theorems and proofs, all authors contributed to the writing of the paper.

Funding. This work was supported by the Engineering and Physical Sciences Research Council grant EP/R020205/1 to Alain Goriely.

Appendix A Cell diagrams and electrical analogy in the mathematical puzzle of Squaring the Square

A similar concept to cell diagrams was devised in the 1940s by four undergraduate students at Trinity College Cambridge, who were trying to find perfect squared squares. The idea is to dissect a square into smaller squares which do not repeat in size. Prior to the work of R.L. Brooks, C.A.B. Smith, A.H. Stone and W.T. Tutte which was published in the classic paper The dissection of rectangles into squares Brooks1940 , a number of squared rectangles and imperfect squares (requiring repeated triangles) had already been known. The four students searched for perfect squared squares mostly empirically, but developed an extensive theory of squared rectangles which combined the theory of planar graphs and of electrical networks, and to which which cell diagrams and cell graphs have a lot of similarity. Their networks, which they called polar nets or p-nets, but which were later dubbed Smith diagrams after C.A.B. Smith, were simply connected polar directed acyclic graphs same as the cell diagrams presented in this work. The students took advantage of Kirchhoff’s laws purely as a tool to balance lengths, interpreting the sides of squares as currents, and potential differences as lengths (with the assumption of unit resistance, and disregarding physical units). Smith diagrams were a semi-empirical tool to organise the imperfect squared squares found in the process towards discovering perfect squared squares, as described later in an essay by W.T. Tutte Tutte1961 . Later an exhaustive computer search revealed that the smallest squared square contains 21 squares, and this is the current logo of the Trinity Mathematical Society111Trinity Mathematical Society logo and its history: https://tms.soc.srcf.net/about-the-tms/the-squared-square/#2. Figure 9 shows the first published squared rectangle known as Lady Isabel’s Casket, as well as the original illustration of p-nets in Brooks1940 . For a detailed history of squared squares and an overview over the use of p-nets, see Anderson2013 . For a concise popular overview see Stewart1995 .