Coupling Chemotaxis and Growth Poromechanics for the Modelling of Feather Primordia Patterning

Abstract

In this paper we propose a new mathematical model for describing the complex interplay between skin cell populations with fibroblast growth factor and bone morphogenetic protein, occurring within deformable porous media describing feather primordia patterning. Tissue growth, in turn, modifies the transport of morphogens (described by reaction-diffusion equations) through diverse mechanisms such as advection from the solid velocity generated by mechanical stress, and mass supply. By performing an asymptotic linear stability analysis on the coupled poromechanical-chemotaxis system (assuming rheological properties of the skin cell aggregates that reside in the regime of infinitesimal strains and where the porous structure is fully saturated with interstitial fluid and encoding the coupling mechanisms through active stress) we obtain the conditions on the parameters—especially those encoding coupling mechanisms—under which the system will give rise to spatially heterogeneous solutions. We also extend the mechanical model to the case of incompressible poro-hyperelasticity and include the mechanisms of anisotropic solid growth and feedback by means of standard Lee decompositions of the tensor gradient of deformation. Because the model in question involves the coupling of several nonlinear PDEs, we cannot straightforwardly obtain closed-form solutions. We therefore design a suitable numerical method that employs backward Euler time discretisation, linearisation of the semidiscrete problem through Newton–Raphson’s method, a seven-field finite element formulation for the spatial discretisation, and we also advocate the construction and efficient implementation of tailored robust solvers. We present a few illustrative computational examples in 2D and 3D, briefly discussing different spatio-temporal patterns of growth factors as well as the associated solid response scenario depending on the specific poromechanical regime. Our findings confirm the theoretically predicted behaviour of spatio-temporal patterns, and the produced results reveal a qualitative agreement with respect to the expected experimental behaviour. We stress that the present study provides insight on several biomechanical properties of primordia patterning.

1. Introduction

Transport of species in tissue is a key process ubiquitous in biological systems. In particular, chemotaxis models [1] describe the directed movement of cells in response to chemicals (attractants or repellents), and can predict the formation of clustered structures. This mechanism has been observed in a variety of embryogenesis processes, such as gastrulation [2] and feather development [3]. In the latter, we refer to the buds that then give origin to feathers as primordia, and their origin has captured the attention of many theoretical and experimental studies. Mathematical models of different complexity have been used for primordia pattern formation, mainly based on the well-known Keller–Segel and Patlak–Keller–Segel-type chemotaxis equations considering fibroblast growth factor (FGF) and bone morphogenic protein (BMP) [4,5,6].

Poroelasticity refers to mixture models in which a solid phase coexists with (at least) one fluid phase/fluid particles [7]. In classical poroelastic structures composed by solid and fluid, the stress of the mixture is the sum of the interstitial pressure of the fluid and the stress of the (hyper)elastic solid phase, and both the fluid and solid phases are assumed to be incompressible. Biological organs and tissues are naturally porous at the tissue level, as they are heterogeneously composed, for example, of both muscle and blood, or of interstitial cerebrospinal fluid and cells, or collagen fibres and extracellular matrix in systems such as cartilage tissue. This phase separation persists even up to the cellular level, as there is the cytoskeleton and the cytoplasm. For this reason, poroelastic models have become very pervasive in the modelling of soft living tissue, such as in oedema formation [8,9], cardiac perfusion [10,11], lung characterisation [12], and brain injury [13].

In the context of biologically-oriented problems, experiments have shown that the rheology of cytoplasm within living cells exhibits a poroelastic behaviour [14], and in turn, the composition of cells and the extracellular matrix constitutes an overall poromechanical system. Using only fluid models lacks in general the possibility of producing heterogeneous proliferation patterns thanks to the stress localisation shear properties observed in deformable solids or solid-fluid mixtures. The presence of chemical solutes locally modifies morphoelastic properties and these processes can be homogenised to obtain macroscopic models of poroelasticity coupled with advection-reaction-diffusion equations (see e.g., [15,16]). In the proposed model we assume that the poroelastic growth tensor incorporates the contribution from the elongation of other structures in the skin such as collagen fibres. In this work we complement and extend the equations for feather primordia spatial organisation proposed by Painter et al. [6].

With this biological basis, the scope of our work is twofold: on one hand, we extend the existing pattern formation models for primordia, by considering their interaction with the intracellular space in the outset of growth. For this, in the regime of infinitesimal strains we incorporate an active stress in the solid phase which is regulated by morphogens. In the finite-strain regime we use a growth model which is consistent with the remodelling theory from, e.g., [17,18,19]. On the other hand, we perform a thorough stability analysis to understand the coupling mechanisms in the model and to investigate the conditions that give rise to pattern formation. The governing equations are of nonlinear nature, and they involve many vector and scalar fields. Attempting to solve them or to extract useful insight from their structure without further simplifications is very complicated. We first simplify and linearise the coupled system and carry out a thorough stability analysis based on non-trivial homogeneous steady states. We also derive appropriate base-case dispersion relations that serve to understand the coupling mechanisms in the proposed model, and they allow us to investigate the conditions that give rise to pattern formation. For this we follow the classical approach from, e.g., [13,20,21,22], where the so-called Routh–Hurwitz assumptions for stability are verified with the aim of further characterising relevant parametric spaces. These investigations allow us to conclude that an oscillatory Turing instability exists above a critical value of the mechanochemical coupling strength, for example.

It is important to mention that our paper also focuses on the numerical simulation of coupled chemotaxis-poromechanics models in 2D and 3D, which validate and extend the analytical results from the linear stability analysis. Many robust methods exist already for poroelasticity, see, e.g., [11,23,24,25,26,27] and the references therein. On the other hand, the theoretical numerical analysis of the properties of finite element discretisations for coupled poromechanics and two-species chemotaxis has been addressed in [28] using total pressure, and in [29] using conservative schemes; all these works focusing on the regime of linear poroelasticity. The modification of the analysis to accommodate the finite-strain framework has been recently conducted in [8,30]. Here we extend the latter formulations to include other morphogens and additional nonlinearities, and focus on a stronger two-way coupling required for the case of primordia patterns. The development of efficient preconditioners for chemotaxis has not yet been addressed in literature, whereas nonlinear poromechanics has already been studied [31,32] albeit only in a small deformations regime. In this work we provide efficient and scalable for both chemotaxis and large deformations poroelasticity.

We have structured the remainder of this paper in the following manner. Section 2 describes the coupled model for poro-mechano-chemical interactions, restricting the presentation to the regime of linear poroelasticity. We give an adimensional form of the governing equations, making precise boundary and initial conditions. In Section 3 we perform a linear stability analysis addressing pattern formation according to Turing instabilities. We separate the discussion in some relevant cases and derive and portray patterning spaces. Section 4 is devoted to extending the model to the case of nonlinear (finite-strain) poroelasticity and material growth, stating also the coupling with chemotaxis in the undeformed configuration. Some numerical examples are given in Section 5 and we close with a summary and discussion of model extensions and perspectives for future work in Section 6.

2. A Coupled Model of Linear Poroelasticity and Chemotaxis

Let us consider a piece of soft material as a porous medium in ${\mathbb{R}}^d$, $d=2,3$, composed by a mixture of incompressible grains and interstitial fluid, whose description can be placed in the context of the classical Biot consolidation problem (see e.g., [33]). In the absence of gravitational forces, of body loads, and of mass sources or sinks, we seek for each time $t\in (0,t_{\mathrm{final}}]$, the displacement of the porous skeleton, $u\left(t\right):\Omega \to {\mathbb{R}}^d$, and the pore pressure of the fluid, $p\left(t\right):\Omega \to \mathbb{R}$, such that

$$\begin{array}{cc}{\partial }_t\left(C_{0}p+{\alpha }_{BW}\mathrm{div}u\right)-\frac{1}{\eta }\mathrm{div}\{\kappa \nabla p\}=0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill \mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\end{array}$$

(1a)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\sigma ={\sigma }_{\mathrm{poroelast}}+{\sigma }_{\mathrm{act}}\hspace{1em}\hspace{1em}& \mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\rho {\partial }_{tt}\mathit{u}-\mathbf{div}\sigma =\mathbf{0}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hfill & \hfill \mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\end{array}$$

(1c)

where $\kappa \left(x\right)$ is the hydraulic conductivity of the porous medium, $\rho$ is the density of the solid material, $\eta$ is the constant viscosity of the interstitial fluid, $C_0$ is the constrained specific storage coefficient, ${\alpha }_{BW}$ is the Biot–Willis consolidation parameter. In (1b) we are supposing that the poromechanical deformations are also actively influenced by microscopic tension generation. A very simple description is given in terms of active stresses: we assume that the total Cauchy stress contains a passive and an active component, where

$${\sigma }_{\mathrm{poroelast}}=\lambda \left(\mathrm{div}u\right)\mathbf{I}+2\mu \epsilon \left(u\right)-{\alpha }_{BW}p\mathbf{I},$$

(2)

and ${\sigma }_{\mathrm{act}}$ is specified in (5), below. The tensor $\epsilon \left(u\right)=\frac{1}{2}(\nabla u+\nabla u^{⊺})$ is that of infinitesimal strains, $\mathbf{I}$ denotes the second-order identity tensor, and $\mu ,\lambda$ are the Lamé constants (shear and dilation moduli) of the solid structure. Equations (1a)–(1c) represent the conservation of mass, the constitutive relation, and the conservation of linear momentum, respectively.

In addition, let us consider a modified Patlak–Keller–Segel model for the distribution of chemotactic cell populations of mesenchymal cells, m, epithelium activation state, e, fibroblast growth factor (FGF), f, and bone morphogenetic protein (BMP), b. The base-line model has been developed in [6] and (after being properly modified to account for the motion of the underlying deformable porous media) it can be summarised as follows

$$\begin{array}{cc}{\partial }_{t}m+{\partial }_{t}u·\nabla m-\mathrm{div}\left(D_m\nabla m-\alpha mexp(-\gamma m)\nabla f\right)=0\hfill & \hfill \hspace{1em}\hspace{1em}\mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\end{array}$$

(3a)

$$\begin{array}{cc}{\partial }_{t}e-[{\kappa }_{1}w(x,t)h_1\left(m\right)+{\kappa }_2{h}_2\left(m\right)](1-e)+[1-h_1\left(m\right)]({\kappa }_3+{\kappa }_{4}b)e=0\hfill & \hfill \hspace{1em}\hspace{1em}\mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\end{array}$$

(3b)

$$\begin{array}{cc}\hfill {\partial }_{t}f+{\partial }_{t}u·\nabla f-\mathrm{div}(D_f\nabla f)-{\kappa }_{F}e+{\delta }_{F}f+{\xi }_f\mathrm{div}u=0& \hfill \hspace{1em}\hspace{1em}\mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\end{array}$$

(3c)

$$\begin{array}{cc}\hfill {\partial }_{t}b+{\partial }_{t}u·\nabla b-\mathrm{div}(D_b\nabla b)-{\kappa }_B{h}_3\left(m\right)m+{\delta }_{B}b=0& \hfill \hspace{1em}\hspace{1em}\mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\end{array}$$

(3d)

where $D_m,D_f,D_b$ are positive definite diffusion matrices, and the spatio-temporal and nonlinear coefficients (in this case, the priming wave responsible for the generation of spatial patterns and the degree of clustering of mesenchymal cells, respectively) are defined as

$$w(x,t)=\frac{{\omega }_1}{2}\{1+tanh\left({\omega }_2[t-x_2/{\omega }_3]\right)\},h_i\left(m\right)=m^{P_i}{[K_i^{P_i}+m^{P_i}]}^{-1},$$

(4)

where ${\kappa }_i,{\omega }_i,P_i,\gamma ,{\delta }_i,{\xi }_f$ are positive model constants. Note that the mechano-chemical feedback (the process where mechanical forces modify the reaction-diffusion effects) is here assumed only through an additional reaction term in the FGF Equation (3c), depending linearly on volume change. Since ${\xi }_f>0$, this contribution acts as a local sink of FGF during dilation.

On the other hand, we assume that the active stress component acts isotropically on the medium (see e.g., [17]), and it depends nonlinearly on the concentration of mesenchymal cells, as proposed for instance in [21]

$${\sigma }_{\mathrm{act}}=\left(\lambda +\frac{2\mu }{3}\right)\frac{\tau m}{1+\zeta m^2}\mathbf{I},$$

(5)

with $\tau$ a model constant to be specified later on, which can be positive (implying that the function modifies the motion of the medium as a local dilation) or negative (isotropically distributed compression).

In order to reduce the number of model parameters, we focus on a dimensionless counterpart of systems (1a)–(1c) and (3a)–(3d), which can be derived using the following transformation, suggested in [6]

$$\begin{array}{c}m=\frac{m^{*}}{\gamma },e=e^{*},f=\frac{{\kappa }_F{f}^{*}}{{\delta }_B},b=\frac{{\kappa }_B{b}^{*}}{\gamma },\\ u=\sqrt{\frac{D_b}{{\delta }_B}}{u}^{*},p=p^{*},t=\frac{t^{*}}{{\delta }_B},x=\sqrt{\frac{D_b}{{\delta }_B}}{x}^{*},\end{array}$$

where hereafter the stars are dropped for notational convenience. The adimensional coupled system is then equipped with appropriate initial data at rest

$$\begin{array}{c}m\left(0\right)=m_0,e\left(0\right)=e_0=\frac{{\kappa }_{1}w(x,0)h_1\left(m_0\right)+{\kappa }_2{h}_2\left(m_0\right)}{{\kappa }_{1}w(x,0)h_1\left(m_0\right)+{\kappa }_2{h}_2\left(m_0\right)+[1-h_1\left(m_0\right)]({\kappa }_3+{\kappa }_{4}b)},\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}f\left(0\right)=\frac{e_0}{{\delta }_F},b\left(0\right)=h_3\left(m_0\right)m_0,u\left(0\right)=\mathbf{0},{\partial }_{t}u\left(0\right)=\mathbf{0},p\left(0\right)=0,\end{array}$$

(6)

defined in $\Omega$; and boundary conditions in the following manner

$$\begin{array}{ccc}\hfill \{D_m\nabla m-\alpha mexp(-\gamma m)\nabla f\}·n=D_f\nabla f·n=\nabla b·n& =0\hfill & \hfill \mathrm{on}\partial \Omega \times (0,t_{\mathrm{final}}],\end{array}$$

(7a)

$$\begin{array}{cc}\hfill u=u_{\mathsf{\Gamma }}\mathrm{and}\frac{\kappa }{\eta }\nabla p·n=0& \hfill \mathrm{on}\mathsf{\Gamma }\times (0,t_{\mathrm{final}}],\end{array}$$

(7b)

$$\begin{array}{ccc}\hfill \sigma =t\mathrm{and}p& =p_0\hfill & \hfill \mathrm{on}\mathsf{\Sigma }\times (0,t_{\mathrm{final}}],\end{array}$$

(7c)

where denotes the outward unit normal on the boundary, and  $\partial \Omega =\mathsf{\Gamma }\cup \mathsf{\Sigma }$ is disjointly split into $\mathsf{\Gamma }$ and $\mathsf{\Sigma }$ where we prescribe clamped boundaries and zero fluid normal fluxes; and a given traction $t$ together with constant fluid pressure $p_0$, respectively.

3. Linear Stability Analysis and Dispersion Relation

3.1. Preliminaries

We proceed to derive a linear stability analysis for the coupled problem (1a)–(5). As usual the analysis is performed on an infinite domain in ${\mathbb{R}}^d$, with $d=2,3$. The first step consists in linearising the poro-mechano-chemical system around a steady state, defined in (6). The linearised dimensionless equations are given by

$$\begin{array}{cc}\hfill {\partial }_t\left(C_{0}p+{\alpha }_{BW}\mathrm{div}u\right)-\frac{1}{\eta }\mathrm{div}\{\kappa \nabla p\}& =0,\hfill \\ \hfill \rho {\partial }_{tt}u-\mathbf{div}\left({\sigma }_{\mathrm{poroelast}}+{\sigma }_{\mathrm{act}}^{\mathrm{lin}}\right)& =\mathbf{0},\hfill \\ \hfill {\partial }_{t}m-\mathrm{div}\{D_m\nabla m-\alpha m_0{\mathrm{e}}^{-m_0}\nabla f\}& =0,\hfill \\ \hfill {\partial }_{t}e-A_m(m_0,e_0)m+A_e(m_0,b_0)e+A_b(m_0,e_0)b& =0,\hfill \\ \hfill {\partial }_{t}f-\mathrm{div}\{D_f\nabla f\}-e+{\delta }_{F}f+{\xi }_f\mathrm{div}u& =0,\hfill \\ \hfill {\partial }_{t}b-\Delta b-H_3\left(m_0\right)m+{\delta }_{B}b& =0,\hfill \end{array}$$

(8)

where

$${\sigma }_{\mathrm{act}}^{\mathrm{lin}}=\left(\left(\lambda +\frac{2\mu }{3}\right)\tau \frac{1-\zeta m_0^2}{{\left(1+\zeta m_0^2\right)}^2}\right)m\mathbf{I},$$

and $A_m$, $A_e$, $A_b$, $H_3$ are functions of the steady state values, which will be made precise in Definition 1, below. We look then for solutions of the form $u,p,m,e,f,b\propto e^{\mathrm{i}k·x+\varphi t}$, where $k$ is the wave vector (a measure of spatial structure) and $\varphi$ is the linear growth factor. By substituting this ansatz on the system (8), we get a system of linear equations for the vector $w={\left(u,p,m,e,f,b\right)}^{⊺}$, where the associated complex eigenvalues $\varphi$, give information on occurrence of instability of the steady state (i.e., pattern formation), when its real component is positive. It is worth pointing out that the present analysis will only address Turing patterning (Hopf bifurcations or other forms of instability that relate to analysing the imaginary part of $\varphi$, are not considered). In order to specify the aforementioned system of equations for $w$, we collect in Proposition 1 some useful preliminary relations. The proof is postponed to the Appendix A.

 Proposition 1. 

Let g, v be sufficiently regular scalar and vector functions defined by $g=exp(\mathrm{i}k·x+\varphi t)$ and $v=v_{0}exp(\mathrm{i}k·x+\varphi t)$ (with $v_0$, a constant vector), respectively. Let us also set $\theta =\mathrm{div}v$. Then

$$\nabla f=\mathrm{i}kf,\Delta f={\mathrm{i}}^{2}fk·k=-k^{2}f,{\partial }_{t}f=\varphi f,\nabla v=\mathrm{i}v\otimes k,\mathrm{div}v=\mathrm{i}v·k,{\partial }_{t}v=\varphi v,$$

(9a)

$$\mathbf{div}\epsilon \left(v\right)=-\frac{(v·k)k+k^{2}v}{2},{\partial }_t\epsilon \left(v\right)=\varphi \epsilon \left(v\right),\mathbf{div}\left(\theta \mathbf{I}\right)=-(v·k)k,{\partial }_t\theta =\varphi \theta .$$

(9b)

The sought system is defined next, starting from (8).

 Definition 1. 

Let $w={\left(u,p,m,e,f,b\right)}^{⊺}\in {\mathbb{R}}^{d+5}$, $d=2,3$ be the vector of independent variables. Then the associated linear system is given by $\mathbf{M}w={\mathbf{0}}_{d+5}$, where the matrix ${\mathbf{M}}_{(d+5)\times (d+5)}$ adopts the form

$$\mathbf{M}=\left[\begin{array}{cc}{\mathbf{M}}_{11}& {\mathbf{M}}_{12}\\ {\mathbf{M}}_{21}& {\mathbf{M}}_{22}\end{array}\right],$$

with the blocks defined as

$$\begin{array}{cc}\hfill {\mathbf{M}}_{11}& =\left[\begin{array}{ccccc}{A}_1{k}_1+B& A_1{k}_2& \cdots & A_1{k}_d& \mathrm{i}{k}_1\\ A_2{k}_1& A_2{k}_2+B& \cdots & A_2{k}_d& \mathrm{i}{k}_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ A_d{k}_1& A_1{k}_2& \cdots & A_d{k}_d+B& \mathrm{i}{k}_d\\ \mathrm{i}{\alpha }_{BW}\varphi k_1& \mathrm{i}{\alpha }_{BW}\varphi k_2& \cdots & \mathrm{i}{\alpha }_{BW}\varphi k_d& C\end{array}\right],\hfill \\ \hfill {\mathbf{M}}_{12}& =\left[\begin{array}{cccc}-\mathrm{i}\left(\lambda +\frac{2\mu }{3}\right)\tau \frac{1-\zeta m_0^2}{{\left(1+\zeta m_0^2\right)}^2}{k}_1& 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮\\ -\mathrm{i}\left(\lambda +\frac{2\mu }{3}\right)\tau \frac{1-\zeta m_0^2}{{\left(1+\zeta m_0^2\right)}^2}{k}_d& 0& 0& 0\end{array}\right],\hfill \\ \hfill {\mathbf{M}}_{21}& =\left[\begin{array}{ccc}0& \cdots & 0\\ 0& \cdots & 0\\ \mathrm{i}{\xi }_f{k}_1& \cdots & \mathrm{i}{\xi }_f{k}_d\\ 0& \cdots & 0\end{array}\right],\hfill \\ \hfill {\mathbf{M}}_{22}& =\left[\begin{array}{cccc}\varphi +D_m{k}^2& 0& -\alpha m_0{\mathrm{e}}^{-m_0}{k}^2& 0\\ -A_m& \varphi +A_e& 0& A_b\\ 0& -1& \varphi +D_f{k}^2+{\delta }_F& 0\\ -H_3& 0& 0& \varphi +k^2+1\end{array}\right],\hfill \end{array}$$

and where the relevant entries are defined as $B=\rho {\varphi }^2+\mu k^2$, $C=C_0\varphi +\frac{\kappa }{\eta }{k}^2$, and 

$$\begin{array}{cc}\hfill A_i& =\left(\mu +\lambda \right)k_i,\mathrm{for}\mathrm{all}i\mathrm{in}\{1,\dots ,d\},\hfill \\ \hfill A_m& =\left({\kappa }_1{\omega }_1{h}_1^{\prime }\left(m_0\right)+{\kappa }_2{h}_2^{\prime }\left(m_0\right)\right)(1-e_0)+h_1^{\prime }\left(m_0\right)({\kappa }_3+{\kappa }_4{b}_0)e_0,\hfill \\ \hfill A_e& ={\kappa }_1{\omega }_1{h}_1\left(m_0\right)+{\kappa }_2{h}_2\left(m_0\right)+(1-h_1\left(m_0\right))({\kappa }_3+{\kappa }_4{b}_0),\hfill \\ \hfill A_b& =(1-h_1\left(m_0\right)){\kappa }_4{e}_0,H_3=h_3\left(m_0\right)+h_3^{\prime }\left(m_0\right)m_0.\hfill \end{array}$$

We then proceed to obtain a dispersion relation associated with the characteristic polynomial $P\left(\varphi \right)=det\left(\mathbf{M}\right)$ of the matrix described in Definition 1. We obtain

$$P(\varphi ;k^2)=B{(\varphi ;k^2)}^{d-1}\left(P_1(\varphi ;k^2)+P_2(\varphi ;k^2)P_3(\varphi ;k^2)\right),$$

(10)

where $B(\varphi ;k^2)=\rho {\varphi }^2+\mu k^2$ is a polynomial with pure imaginary roots. Consequently, it does not have an influence on the stability of the steady state of system (1a)–(5). The polynomials $P_i$, $i=1,2,3$, are given in what follows.

 Definition 2. 

The polynomials conforming the characteristic Equation (10) are

$$\begin{array}{cc}\hfill P_1(\varphi ;k^2)& =b_3{\varphi }^3+b_2{\varphi }^2+b_1\varphi +b_0,\hfill \\ \hfill P_2(\varphi ;k^2)& =a_4{\varphi }^4+a_3{\varphi }^3+a_2{\varphi }^2+a_1\varphi +a_0,\hfill \\ \hfill P_3(\varphi ;k^2)& =c_3{\varphi }^3+c_2{\varphi }^2+c_1\varphi +c_0,\hfill \end{array}$$

(11)

where all coefficients in (11) adopt the following forms

$$\begin{array}{cc}\hfill a_0& =D_m{D}_f(k_{\mathrm{on}}+k_{\mathrm{off}})k^6+\left((D_m{\delta }_B+D_m{D}_f)(k_{\mathrm{on}}+k_{\mathrm{off}})-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m\right)k^4\hfill \\ \hfill & +\left(D_m{\delta }_B(k_{\mathrm{on}}+k_{\mathrm{off}})+\alpha m_0{\mathrm{e}}^{-m_0}{H}_3{A}_b-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m\right)k^2,\hfill \\ \hfill a_1& =D_m{D}_f{k}^6+\left((D_m{D}_f+D_m+D_f)(k_{\mathrm{on}}+k_{\mathrm{off}})+D_m{\delta }_B+D_{m}Df\right)k^4\hfill \\ \hfill & +\left((D_m{\delta }_B+{\delta }_B+D_m+D_f)(k_{\mathrm{on}}+k_{\mathrm{off}})+D_m{\delta }_B-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m\right)k^2+{\delta }_B(k_{\mathrm{on}}+k_{\mathrm{off}}),\hfill \\ \hfill a_2& =\left(D_m{D}_f+D_m+D_f\right)k^4+\left(D_m{\delta }_B+(D_m+D_f+1)(k_{\mathrm{on}}+k_{\mathrm{off}})+\delta +D_m+D_f\right)k^2\hfill \\ \hfill & +({\delta }_B+1)(k_{\mathrm{on}}+k_{\mathrm{off}})+{\delta }_B,\hfill \\ \hfill a_3& =\left(D_m+D_f+1\right)k^2+{\delta }_B+k_{\mathrm{on}}+k_{\mathrm{off}}+1,a_4=1,\hfill \\ \hfill b_i& =-\left({\xi }_f\left(\lambda +\frac{2\mu }{3}\right)\tau \frac{1-\zeta m_0^2}{{\left(1+\zeta m_0^2\right)}^2}\alpha m_0{\mathrm{e}}^{-m_0}{k}^2\right){\widehat{b}}_i,{\widehat{b}}_0=\frac{\kappa }{\eta }\left(k_{\mathrm{on}}+k_{\mathrm{off}}\right)\left(k^2+1\right)k^2,\hfill \\ \hfill {\widehat{b}}_1& =\frac{\kappa }{\eta }{k}^4+\left(\left(C_0+\frac{\kappa }{\eta }\right)\left(k_{\mathrm{on}}+k_{\mathrm{off}}\right)+\frac{\kappa }{\eta }\right)k^2+C_0\left(k_{\mathrm{on}}+k_{\mathrm{off}}\right),\hfill \\ \hfill {\widehat{b}}_2& =\left(C_0+\frac{\kappa }{\eta }\right)k^2+C_0\left(k_{\mathrm{on}}+k_{\mathrm{off}}+1\right),{\widehat{b}}_3=C_0,\hfill \\ \hfill c_0& =\frac{\kappa }{\eta }(2\mu +\lambda )k^4,c_1=\left(C_0(2\mu +\lambda )+{\alpha }_{BW}\right)k^2,c_2=\rho \frac{\kappa }{\eta }{k}^2,c_3=\rho C_0,\hfill \\ \hfill k_{\mathrm{on}}& ={\kappa }_1{\omega }_1{h}_1\left(m_0\right)+{\kappa }_2{h}_2\left(m_0\right),k_{\mathrm{off}}=(1-h_1\left(m_0\right))({\kappa }_3+{\kappa }_4{b}_0).\hfill \end{array}$$

From Definition 2 we immediately see that the characteristic polynomial is of up to seventh degree, making it difficult to determine analytically the main features of the system.

Next we concentrate on addressing some specific scenarios of particular interest. Unless specified otherwise, we will employ the parameter values provided in Box 1. Note that the term ${\sigma }_{\mathrm{act}}^{\mathrm{lin}}$ vanishes for $m_0=1$, so, differently from [6], we modify the steady state for mesenchymal cells concentration to $m_0=2$.

Box 1. Model parameters used in the linear stability analysis of the poro-mechano-chemical system

$m_0=2,D_m=0.01,D_f=0.1,\alpha =4,{\kappa }_1=0.05,{\kappa }_2=0.025,{\kappa }_3=1,{\kappa }_4=1,$

$K_1=1,K_2=2,K_3=5,P_1=P_2=P_3=2,{\delta }_F=1,{\omega }_1=1,{\omega }_2=5,{\omega }_3=0.04$

$E=3·{10}^4,\nu =0.4,C_0={10}^{-3},\kappa ={10}^{-4},{\alpha }_{BW}=0.1,\tau =60,\eta =0.1,{\xi }_f=0.15,\rho =1,\zeta =1$

3.2. Spatially Homogeneous Distributions

If $k^2=0$ then the characteristic polynomial $P(\varphi ;k^2)$ is

$$\begin{array}{cc}\hfill P(\varphi ;0)=& \rho C_0{\varphi }^8[{\varphi }^3+\left({\delta }_B+k_{\mathrm{on}}+k_{\mathrm{off}}+1\right){\varphi }^2\hfill \\ \hfill & +\left(({\delta }_B+1)(k_{\mathrm{on}}+k_{\mathrm{off}})+{\delta }_B\right)\varphi +\delta (k_{\mathrm{on}}+k_{\mathrm{off}})].\hfill \end{array}$$

The well-known Routh–Hurwitz conditions (see e.g., [34]) state that for any polynomial of order 3, a necessary and sufficient set of conditions should be satisfied in order to ensure that all the roots are in the space of complex non-positive real values, ${\mathbb{C}}_{-}=\{z\in \mathbb{C}:\Re \left(z\right)\le 0\}$. For a general polynomial $P\left(\varphi \right)={\alpha }_3{\varphi }^3+{\alpha }_2{\varphi }^2+{\alpha }_1\varphi +{\alpha }_0$, where all the ${\alpha }_i>0$, we have

$${\alpha }_2{\alpha }_1-{\alpha }_3{\alpha }_0>0.$$

(12)

In our case, ${\alpha }_2{\alpha }_1-{\alpha }_3{\alpha }_0=\left(({\delta }_B+1)(k_{\mathrm{on}}+k_{\mathrm{off}})+{\delta }_B\right)\left({\delta }_B+k_{\mathrm{on}}+k_{\mathrm{off}}\right)+{\delta }_B+k_{\mathrm{on}}+k_{\mathrm{off}}$, which is a real positive constant. This indicates that the steady state is stable in a spatially homogeneous system, irrespective of the value of the model parameters.

3.3. Zero Chemotaxis

When $\alpha =0$, all the coefficients $b_i$ from Definition 2 are zero. Consequently, the characteristic polynomial is defined by the sub-polynomials $P_2$ and $P_3$ only. As all the parameters implied in $c_i$ are positive and the Routh–Hurwitz condition (12) is satisfied, we know that the $P_3$ polynomial has only non-positive real part complex solutions and therefore does not influence the stability of the steady state. Thus, the instability condition can only come from $P_2$, related only to the chemotaxis model analysed in [6]. Substituting $\alpha =0$ in the sub-matrix ${\mathbf{M}}_{22}$ leads to a dispersion relation

$$P\left(\varphi \right)=(\varphi +D_m{k}^2)(\varphi +A_e)(\varphi +D_f{k}^2+{\delta }_F)(\varphi +k^2+1),$$

and consequently the roots are negative real numbers. The absence of chemotaxis prevents then any pattern formation irrespective of the presence or absence of poro-mechanical coupling.

3.4. Uncoupled System

Similarly as in the no-chemotaxis scenario above, imposing $\tau =0$ or ${\xi }_f=0$ leads to $P_1(\varphi ;k^2)=0$ for all $\varphi ,k^2$. Additionally, as $P_3(\varphi ;k^2)$ has only non-positive roots, only the chemotaxis sub-polynomial enables us to find eigenvalues that lead to unstable systems. For $P_2$ the Routh–Hurwitz conditions are given by

$$a_i>0,a_2{a}_3-a_1{a}_4>0\mathrm{and}{a}_1{a}_2{a}_3-a_0{a}_3^2-a_1^2{a}_4>0.$$

(13)

Then, from Definition 2, we can readily observe that $a_4$, $a_3$ and $a_2$ are real positive coefficients for any positive parameter value. A complete analysis of (13) is analytically quite involved. Nevertheless, and in accordance with the analysis in [6], we can restrict the discussion to the conditions that break $a_0>0$.

We first observe, from its definition, that $a_0$ is a polynomial of even order with respect to $k^2$, guaranteeing that it has at least one local extrema. Consequently, we look for the critical wave number, $k_c^2>0$, defined by the equation $a_0^{\prime }\left(k^2\right)=0$, that substituting in $a_0$ will lead exactly to be at zero for the associated critical parameter ${\theta }_c$ that we want to analyse (i.e., $a_0({\theta }_c;k_c^2)=0$). In the uncoupled scenario, $a_0^{\prime }$ is a quadratic polynomial with respect to $k^2$ and roots can be easily obtained. Nonetheless, in order to have $a_0$ negative for positive $k^2$, one has to satisfy the following three conditions

$$\begin{array}{cc}\hfill (D_m{\delta }_B+D_m{D}_f)(k_{\mathrm{on}}+k_{\mathrm{off}})-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m& <0,\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill D_m{\delta }_B(k_{\mathrm{on}}+k_{\mathrm{off}})+\alpha m_0{\mathrm{e}}^{-m_0}{H}_3{A}_b-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m& <0,\hfill \\ \hfill 4{\left((D_m{\delta }_B+D_m{D}_f)(k_{\mathrm{on}}+k_{\mathrm{off}})-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m\right)}^2-& \hfill \end{array}$$

(14b)

$$\begin{array}{cc}\hfill 12D_m{D}_f(k_{\mathrm{on}}+k_{\mathrm{off}})\left(D_m{\delta }_B(k_{\mathrm{on}}+k_{\mathrm{off}})+\alpha m_0{\mathrm{e}}^{-m_0}{H}_3{A}_b-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m\right)& >0.\hfill \end{array}$$

(14c)

Either (14a) or (14b) enforces that $k^2$ has a real positive part, while (14c) is necessary to have a real wave number. These conditions are similar to the ones stated in [6]. The inequality (14b) emphasises how the positive feedback mechanism of the accumulation of mesenchymal cells through chemotaxis leads to stimulate further FGF secretion, based on the activation of the epithelium. When the condition is satisfied, patterning is achieved.

Condition (14a) indicates that the steady state is destabilised when $A_m>0$, occurring when the clustering-mediated activation is sufficiently strong. Figure 1 (panels (A1)–(B1)) present the patterning space based on the implicit functions defined in (14a)–(14c) for the $(m_0,\alpha )$ and $({\kappa }_1,{\kappa }_4)$ space. The conditions are represented in both Figures by a plain-green curve (14a), a red-dot-dashed curve (14b), and a blue-dashed curve (14c). The boundaries present very similar results to the linear analysis from [6]. Here, we accentuate the region by filling it with light grey. Panels (A2)–(A3) and (B2)–(B3) in Figure 1 portray the functions $a_0\left(k^2\right)$ (parameter condition) and $\lambda \left(k^2\right)$ (dispersion relation) for the studied parameters $m_0$ and ${\kappa }_4$. We present in different colours the behaviour of the system for the critical parameter value (green curve), as well as parameter values obtained by increasing and decreasing by 25% and 50% the critical value, where we recall that fixed parameters are taken from Box 1. We see that moving above or below the critical value results in a finite interval of wave numbers for which the complex eigenvalue presented positive real component and thus leading to instability (patterning). As illustrated also in [6], for a specific parameter set, panels (A1)–(A3) from Figure 1 show that a sufficiently dense dermis cells concentration is necessary to produce instability. While increasing the chemotaxis sensibility of the system we spread the range of possible mesenchymal cells density, by the u-shape form of the boundary, a sufficiently high density of cells can prevent pattern formation (see Figure 1(A1)). The capacity of BMP to deactivate the epithelium by increasing ${\kappa }_4$, rapidly decreases the patterning ability of the system.

3.5. Zero Activation/Inactivation of Epithelium

Different scenarios are possible in such case: only zero activation (${\kappa }_1={\kappa }_2=0$, $e_0=0$), only zero inactivation (${\kappa }_3={\kappa }_4=0$, $e_0=1$) or both zero activation and inactivation of the epithelium. Contrary to [6], the coupling with the poroelastic structure leads to a more involved analysis of the system and, in particular, it does not guarantee emerging patterns. In any case, the coefficients $A_m$ and $A_b$ will be zero and bring the coefficients of $P_2$ to be strictly positive. As the coefficients of $P_3$ are also non-negative, $P(\varphi ;k^2)$ can only lead to negative eigenvalues if the coefficients of $P_1$ are negative. By Definition 2, we know that the $b_i$’s are negative only if the coupled-dependent constant is itself negative, and this phenomenon occurs if one of the following conditions is satisfied

$$-\left(\lambda +\frac{2\mu }{3}\right)\tau \frac{1-\zeta m_0^2}{{\left(1+\zeta m_0^2\right)}^2}>0\iff \{\begin{array}{ccc}\tau >0& \mathrm{and}& m_0<\sqrt{\frac{1}{\zeta }}\\ \tau <0& \mathrm{and}& m_0>\sqrt{\frac{1}{\zeta }}\end{array}.$$

If one of the condition is satisfied, the formation of patterns is expected. This is an opposite conclusion as that drawn in the case of chemotaxis only [6].

3.6. General Case

The thorough analysis of the coupled system involves a seventh-order polynomial. Analytical Routh–Hurwitz conditions in this case are hardly accessible and the crossed effects due to multiple parameters and conditions will inevitably make difficult the interpretation of the system behaviour. We can however restrict then the analysis to ${\varphi }^0$, with coefficients $d_0=a_0{c}_0+b_0$, where $a_0$, $b_0$, $c_0$ are as in Definition 2. After some algebraic manipulations, we derive $d_0$ with respect to $k^2$ and obtain

$$d_0^{\prime }\left(k^2\right)={\theta }_4{k}^8+{\theta }_3{k}^6+{\theta }_2{k}^4+{\theta }_1{k}^2,$$

where

$$\begin{array}{cc}\hfill {\theta }_1& =-2{\xi }_f\left(\lambda +\frac{2\mu }{3}\right)\tau \frac{1-\zeta m_0^2}{{\left(1+\zeta m_0^2\right)}^2}\alpha m_0{\mathrm{e}}^{-m_0}(k_{\mathrm{on}}+k_{\mathrm{off}}),\hfill \\ \hfill {\theta }_2& =3[(2\mu +\lambda )\left(D_m{\delta }_B(k_{\mathrm{on}}+k_{\mathrm{off}})+\alpha m_0{\mathrm{e}}^{-m_0}{H}_3{A}_b-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m\right)\hfill \\ \hfill & -{\xi }_f\left(\lambda +\frac{2\mu }{3}\right)\tau \frac{1-\zeta m_0^2}{{\left(1+\zeta m_0^2\right)}^2}\alpha m_0{\mathrm{e}}^{-m_0}(k_{\mathrm{on}}+k_{\mathrm{off}})],\hfill \\ \hfill {\theta }_3& =4(2\mu +\lambda )\left((D_m{\delta }_B+D_m{D}_f)(k_{\mathrm{on}}+k_{\mathrm{off}})-\alpha m_0{\mathrm{e}}^{-m_0}{A}_m\right),\hfill \\ \hfill {\theta }_4& =5D_m{D}_f(2\mu +\lambda )(k_{\mathrm{on}}+k_{\mathrm{off}}).\hfill \end{array}$$

In contrast with the uncoupled case studied above, here $k_c^2$ is found by solving a quartic polynomial (employing the intrinsic MATLAB function roots, after constructing the associated vector of polynomial coefficients). In addition, one needs to satisfy some additional conditions leading to $d_0<0$ for any positive $k^2$. More specifically, the constraints consist of at least one of the Routh–Hurwitz conditions [34], stated as

$$\begin{array}{cc}\hfill {\theta }_1>0,{\theta }_2>0,{\theta }_3& >0,\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill {\theta }_2{\theta }_3-{\theta }_1{\theta }_4& >0,\hfill \end{array}$$

(15b)

$$\begin{array}{cc}\hfill {\theta }_1{\theta }_2{\theta }_3-{\theta }_1^2{\theta }_4& >0,\hfill \end{array}$$

(15c)

or the discriminant

$$-27{\theta }_4^2{\theta }_1^4+18{\theta }_4{\theta }_3{\theta }_2{\theta }_1^3-4{\theta }_4{\theta }_2^3{\theta }_1^3-4{\theta }_3^3{\theta }_1^3+{\theta }_3^2{\theta }_2^2{\theta }_1^2>0.$$

(16)

Condition (15b) highlights the effect of poromechanics towards the instability of the system, in comparison against the uncoupled cased. Increasing the influence of the mesenchymal cells through $\tau$, or the influence of mechanics on f through ${\xi }_f$, can relax the requirement that the feedback-positive loop from (14b) should be sufficiently strong to create instability. As mentioned in the activation/inactivation of epithelium scenario, this will depend mostly on the sign of the active stress component of the coefficient, and thus on the sign of $\tau$ and the value of $m_0$, see Figure 2(C1). The dispersion relation for the $(m_0,\alpha )$- and $({\kappa }_1,{\kappa }_2)$– parameter spaces illustrates how, for the same fixed set of parameters from Box 1, the poromechanical coupling leads to a decrease of the critical parameter value (Figure 2(A3,B3)). Unfortunately, the ease to push the system to instability is counter-balanced by a more restrictive space of parameters values where patterns can be formed (more precisely, compared to Figure 1(A1,B1) and Figure 2(A1,B1,C1)). This can increase the difficulty in tuning parameters to produce specific patterns. Indeed, zooming into the parameter pair $(m_0,\alpha )$ on the region where the boundary conditions cross (see Figure 3), exposes how easily the patterns can disappear after a relatively small variation in parameter values.

4. Extension to Finite-Strain Poroelasticity and Growth

We briefly recall some kinematic considerations, needed in order to incorporate growth effects. Thorough reviews of mechanical models for growth of soft living tissues can be found in [17,18]. Let us denote by $\mathbf{X}$ the position of a material point in the reference (or initial) configuration (region ${\Omega }_0$ with boundary $\partial {\Omega }_0={\mathsf{\Gamma }}_0\cup {\mathsf{\Sigma }}_0$ and outward pointing unit normal ${\mathbf{N}}_0$) and let $x=\phi \left(\mathbf{X}\right)$ denote its position at time t in the current configuration of the domain, $\Omega$, determined by the deformation mapping $\phi (·,t):{\Omega }_0\to \Omega$. We adopt the notation $\mathbf{Div},\mathrm{Div}$ for the divergence of tensor and vector fields, and $\mathbf{Grad},\mathrm{Grad}$ for the gradients of vector and scalar fields, respectively; all with respect to the material coordinates. The geometric deformation tensor is $\mathbf{F}=\mathbf{Grad}\phi =\mathbf{I}+\mathbf{Grad}u$ and its Jacobian is $J=det\mathbf{F}>0$. The velocity in Lagrangian coordinates is denoted as $v=\dot{u}$, and $\mathbf{C}={\mathbf{F}}^{⊺}\mathbf{F}$ is the right Cauchy–Green strain tensor.

As in the linearised regime, we assume that the region occupied by the tissue is fully saturated, and consisting of a fluid and of a solid phase:

$${\varphi }_f+{\varphi }_s=1,$$

(17)

where the ${\varphi }_j$’s denote the volume fraction of each phase. The phase s is assumed a neo–Hookean solid whereas the liquid phase is considered as an ideal fluid. Both phases are assumed intrinsically incompressible. The expected growth is due to proliferation of cells, and implying a mass exchange between the phases. The mass balances in the Eulerian framework assume the form

$$\begin{array}{cc}\hfill \frac{\partial {\varphi }_s}{\partial t}+\mathrm{div}\left({\varphi }_s{v}_s\right)& =r_s\mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\hfill \end{array}$$

(18a)

$$\begin{array}{cc}\hfill \frac{\partial {\varphi }_f}{\partial t}+\mathrm{div}\left({\varphi }_f{v}_f\right)& =-r_s\mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\hfill \end{array}$$

(18b)

with $v_s,v_f$ the Eulerian velocities in each phase, and $r_s={\alpha }_0{\varphi }_s(1-{\varphi }_s)(m-m_0)$ being the mass growth rate (net proliferation rate of cells in the solid phase), here assumed linearly dependent on the concentration of mesynchemal cells and decreasing as the available space decreases. In turn, the general form of the momentum balances is

$$\begin{array}{cc}\hfill \mathbf{div}\left({\varphi }_s{\sigma }_{\mathrm{eff},s}\right)-{\varphi }_s\nabla p+f_{sf}& =\mathbf{0}\mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\hfill \\ \hfill \mathbf{div}\left({\varphi }_f{\sigma }_{\mathrm{eff},f}\right)-{\varphi }_f\nabla p+f_{fs}& =\mathbf{0}\mathrm{in}\Omega \times (0,t_{\mathrm{final}}],\hfill \end{array}$$

where the ${\sigma }_{\mathrm{eff},j}$’s are the respective effective Cauchy stresses, p is the average pressure of the liquid phase, and 

$$\begin{array}{cc}\hfill f_{sf}& =p\nabla {\varphi }_s+\eta /\kappa {\varphi }_s{\varphi }_f(v_f-v_s),\hfill \\ \hfill f_{fs}& =p\nabla {\varphi }_f-\eta /\kappa {\varphi }_f{\varphi }_s(v_f-v_s),\hfill \end{array}$$

denote drag terms and net sources of momentum, where $\kappa ,\eta$ are the permeability and dynamic viscosity, respectively. The effective stress of the liquid phase ${\sigma }_{\mathrm{eff},f}$ is assumed negligible with respect to the interstitial pressure gradient and therefore the mixture stress is ${\sigma }_{\mathrm{eff}}={\sigma }_{\mathrm{eff},s}$. In addition, Darcy’s law gives

$$v_f-v_s=-\frac{\kappa }{\eta {\varphi }_f}\nabla p,$$

(19)

and then the mixture momentum balance is written as

$$-\mathbf{div}\left({\sigma }_{\mathrm{eff}}\right)+\nabla p=\mathbf{0}\mathrm{in}\Omega \times (0,t_{\mathrm{final}}].$$

(20)

The deformation mapping can be decomposed into a purely growth part and the remainder, elastic, deformation. Such splitting implies that there exists an intermediate configuration $\tilde{\Omega }$ between ${\Omega }_0$ and $\Omega$ which is not necessarily compatible (and which we assume is completely stress free, see sketch in Figure 4), and consequently a multiplicative decomposition of the deformation gradient into a growth deformation gradient (describing local generation or removal of material points) and an elastic deformation gradient tensor is admitted (see, for instance, [35,36,37]).

$$\mathbf{F}={\mathbf{F}}_e{\mathbf{F}}_g,$$

(21)

and therefore $J=J_e{J}_g$ with $J_e:=det{\mathbf{F}}_e$, $J_g:=det{\mathbf{F}}_g$. Note also that since $\mathbf{F}$ is assumed nonsingular, so are the tensors ${\mathbf{F}}_e,{\mathbf{F}}_g$.

Since the intermediate configuration is considered stress-free, stresses are exerted only by the elastic deformation, and the constitutive relations between a given strain energy function $\Psi$ (that characterises the material response of the solid for hyperelastic materials) and the measures of stress can be stated with respect to the intermediate configuration $\tilde{\Omega }$. For the Cauchy stress this gives

$${\sigma }_{\mathrm{eff}}=2J_e^{-1}{\mathbf{F}}_e\frac{\partial \Psi }{\partial {\mathbf{C}}_e}{\mathbf{F}}_e^{-⊺},$$

(22)

where ${\mathbf{C}}_e={\mathbf{F}}_e^{⊺}{\mathbf{F}}_e$ and $J_e=det\left(\mathbf{F}{\mathbf{F}}_g^{-1}\right)$ is the reversible part of the Jacobian. This constitutive equation requires to describe also $\mathbf{F}$ and ${\mathbf{F}}_g$, but from the relation $\dot{\mathbf{F}}{\mathbf{F}}^{-1}=\nabla v_s$, we realise that only ${\mathbf{F}}_g$ needs to be specified.

In order to translate Equations (18a)–(22) into Lagrangian form, we recall the change in volume and area relations $\mathrm{d}V=J\mathrm{d}{V}_0$ and $\mathrm{d}S=J{\mathbf{F}}^{-⊺}\mathrm{d}{S}_0$ together with the transformation of unit outward normal vectors $={\mathbf{F}}^{-⊺}{\mathbf{N}}_0$. For the solid mass balance (18a) it suffices to apply Reynolds transport theorem, the change of coordinates, and Maxwell’s localisation to obtain

$$\dot{\overline{J{\varphi }_s}}=J{r}_s.$$

(23)

The Lagrangian form of Darcy’s law

$$v_f-v_s=-\mathbf{F}\frac{\kappa }{\eta {\varphi }_f}{\mathbf{C}}^{-1}\mathrm{Grad}p,$$

(24)

results from (19) in combination with transport, localisation, as well as pull-back operations. For the liquid mass balance (18b), one can apply the generalised Reynolds transport theorem (24), divergence’s theorem, and localisation, to arrive at

$$\dot{\overline{J{\varphi }_f}}+\mathrm{Div}\left(J\frac{\kappa }{\eta }{\mathbf{C}}^{-1}\mathrm{Grad}p\right)=-J{r}_s.$$

(25)

Summing up (23) and (25) and taking an approximation of the fluid volume fraction in terms of fluid pressure and volume change we obtain (28b), below.

For transforming the momentum balance (20) we use the divergence theorem, the change of coordinates, we recall the relation between Cauchy and first Piola–Kirchhoff stress $\mathbf{P}=J\sigma {\mathbf{F}}^{⊺}$, and use again the divergence theorem to get (28a), below. The hyperelastic component to the stress is

$${\mathbf{P}}_e=\frac{\partial \Psi }{\partial \mathbf{F}}{|}_{\tilde{\Omega }},$$

(26)

and for a neo–Hookean material, Equations (21)–(26) imply that the effective stress tensor is

$${\mathbf{P}}_e=J(\mu {\mathbf{B}}_e-\psi \mathbf{I}){\mathbf{F}}^{-⊺},$$

(27)

where $\mu$ is the shear modulus, $\psi$ denotes a Lagrange multiplier, and ${\mathbf{B}}_e={\mathbf{F}}_e{\mathbf{F}}_e^{⊺}$ is the right Cauchy–Green deformation tensor associated with the intermediate configuration $\tilde{\Omega }$. Note that in this context, rather than material incompressibility of the mixture, the Lagrange multiplier $\psi$ enforces a mass conservation kinematic constraint of the form $J_e=1$, confirming that for a hyperelastically incompressible material all volume changes are solely due to growth factors ($J_g=det{\mathbf{F}}_g$ is in turn irreversible) [38]. If $J_g>1$ the body undergoes growth.

In summary, the equations of nonlinear poroelasticity with growth (including momentum conservation of the solid-fluid mixture, fluid mass balance, and incompressibility of elastic deformations) read as follows

$$\begin{array}{ccc}\hfill \rho \ddot{u}-\mathbf{Div}({\mathbf{P}}_e+{\mathbf{P}}_f)& =\mathbf{0}\hfill & \hfill \mathrm{in}{\Omega }_0\times (0,t_{\mathrm{final}}],\end{array}$$

(28a)

$$\begin{array}{ccc}\hfill \dot{\overline{\left(C_{0}p+{\alpha }_{BW}[\mathrm{tr}\mathbf{C}-d]\right)}}-\frac{1}{J}\mathrm{Div}\left(\frac{\kappa }{\eta }J{\mathbf{C}}^{-1}\mathrm{Grad}p\right)& =\ell \hfill & \hfill \mathrm{in}{\Omega }_0\times (0,t_{\mathrm{final}}],\end{array}$$

(28b)

$$\begin{array}{ccc}\hfill J_e& =1\hfill & \hfill \mathrm{in}{\Omega }_0\times (0,t_{\mathrm{final}}],\end{array}$$

(28c)

where the stress exerted by the fluid (when pulled back to the reference configuration) is simply ${\mathbf{P}}_f=-{\alpha }_{BW}Jp{\mathbf{F}}^{-⊺}$, and the divergence and gradient operators are now understood with respect to the undeformed coordinates. In contrast with the model from Section 2, now the fluid has a source that we relate to the growth process and make it precise below. The divergence of displacement accompanying the Biot–Willis constant has now been replaced by $\mathrm{tr}\mathbf{C}-d$, where d is the spatial dimension. The first term on the left-hand side of (28b) represents the rate of change of the total amount of fluid.

It remains to characterise the growth deformation gradient. We suppose that growth occurs due to a set of internal variables ${\gamma }_1,{\gamma }_2,{\gamma }_3$ acting on local orthonormal directions $k_1$, $k_2$, $k_3$ in the undeformed configuration (and in 2D we consider only ${\gamma }_1,{\gamma }_2$ acting on $k_1$, $k_2$, respectively), and a general form for the growth deformation tensor is simply

$${\mathbf{F}}_g=\mathbf{I}+\sqrt{{\gamma }_1}{k}_1\otimes k_1+\sqrt{{\gamma }_2}{k}_2\otimes k_2+{\gamma }_3{k}_3\otimes k_3.$$

Assuming that a transversely isotropic growth due to a constant rate occurs on the plane $(k_1,k_2)$ and the growth due to morphogen concentrations acts mainly on the direction $k_3$, we employ the following specification for the growth factors ${\gamma }_i$:

$${\gamma }_1={\delta }_{1}t,{\gamma }_2={\delta }_{1}t,{\gamma }_3={\delta }_{2}t+{\delta }_3\frac{m}{1+m^2},$$

(29)

where ${\delta }_1,{\delta }_2,{\delta }_3$ are positive constants and we recall that m is the concentration of mesenchymal cells in the form used in the active stress from (5). In the 2D case, (29) will be replaced by

$${\gamma }_1={\gamma }_2={\delta }_{2}t+{\delta }_3\frac{m}{1+m^2}.$$

(30)

As in [19,39] a more general description that defines a functional relation between the rates of the growth factors and the exchange of mass implies a dependence on the solid volume fraction

$$\sum _{i=1}^d\frac{\dot{{\gamma }_i}}{{\gamma }_i}=\frac{r_s}{{\varphi }_s}.$$

Following [30], we can write a constitutive equation for the mass source in (28b) depending on the constants ${\gamma }_i$ as follows

$$\ell ={\ell }_0(2{\delta }_1+{\delta }_2+{\delta }_3{\partial }_m\left({\gamma }_3\right){\partial }_{t}m),$$

(31)

where ${\ell }_0$ is a constant and where we stress that $\rho$ in this section refers to the density of the solid in the reference configuration.

The set of equations is complemented with suitable boundary and initial conditions, formulated in the reference configuration. In contrast with (7b)–(7c), here we impose Robin conditions for the poroelasticity on the whole undeformed boundary

$$\begin{array}{cc}\hfill ({\mathbf{P}}_e+{\mathbf{P}}_f){\mathbf{N}}_0+{\xi }_{R}J{\mathbf{F}}^{-⊺}u& =\mathbf{0}\mathrm{on}({\mathsf{\Gamma }}_0\cup {\mathsf{\Sigma }}_0)\times (0,t_{\mathrm{final}}],\hfill \end{array}$$

(32a)

$$\begin{array}{cc}\hfill \frac{\kappa }{\eta }J{\mathbf{C}}^{-1}\mathrm{Grad}p·{\mathbf{N}}_0& =0\mathrm{on}{\mathsf{\Gamma }}_0\times (0,t_{\mathrm{final}}],p=p_D\mathrm{on}{\mathsf{\Sigma }}_0\times (0,t_{\mathrm{final}}],\hfill \end{array}$$

(32b)

where ${\xi }_R$ is the stiffness of the springs attached to $\partial {\Omega }_0$.

Moreover, the chemotaxis system is recast in the reference domain. In particular, the advection terms are absorbed by the material derivatives (since we are assuming that all species are transported by the whole mixture and can diffuse in both phases) and the diffusion coefficients are modified by the Piola transformation. As in (28b), the divergence of displacement in the sink term from (3c) now is replaced by $\mathrm{tr}\mathbf{C}-d$. The transformation also involves Reynolds transport theorem, divergence theorem, and localisation. Assuming that $D_m,D_f,D_b$ are isotropic, we have

$$\begin{array}{ccc}\hfill \dot{m}-\frac{1}{J}\mathrm{Div}\left(D_{m}J{\mathbf{C}}^{-1}\mathrm{Grad}m-\alpha mexp(-\gamma m)J{\mathbf{C}}^{-1}\mathrm{Grad}f\right)& =0\hfill & \hfill \mathrm{in}{\Omega }_0\times (0,t_{\mathrm{final}}],\end{array}$$

(33a)

$$\begin{array}{ccc}\hfill \dot{e}-[{\kappa }_{1}w(\mathbf{X},t)h_1\left(m\right)+{\kappa }_2{h}_2\left(m\right)](1-e)+[1-h_1\left(m\right)]({\kappa }_3+{\kappa }_{4}b)e& =0\hfill & \hfill \mathrm{in}{\Omega }_0\times (0,t_{\mathrm{final}}],\end{array}$$

(33b)

$$\begin{array}{ccc}\hfill \dot{f}-\frac{1}{J}\mathrm{Div}\left(D_{f}J{\mathbf{C}}^{-1}\mathrm{Grad}f\right)-{\kappa }_{F}e+{\delta }_{F}f+{\xi }_f(\mathrm{tr}\mathbf{C}-d)& =0\hfill & \hfill \mathrm{in}{\Omega }_0\times (0,t_{\mathrm{final}}],\end{array}$$

(33c)

$$\begin{array}{ccc}\hfill \dot{b}-\frac{1}{J}\mathrm{Div}\left(D_{b}J{\mathbf{C}}^{-1}\mathrm{Grad}b\right)-{\kappa }_B{h}_3\left(m\right)m+{\delta }_{B}b& =0\hfill & \hfill \mathrm{in}{\Omega }_0\times (0,t_{\mathrm{final}}],\end{array}$$

(33d)

where divergence and gradients are taken with respect to the undeformed configuration, and the zero-flux boundary conditions (7a) are changed accordingly

$$\begin{array}{cc}\hfill \left(D_{m}J{\mathbf{C}}^{-1}\mathrm{Grad}m-\alpha mexp(-\gamma m)J{\mathbf{C}}^{-1}\mathrm{Grad}f\right)·{\mathbf{N}}_0& =D_{f}J{\mathbf{C}}^{-1}\mathrm{Grad}f·{\mathbf{N}}_0\hfill \\ \hfill & =D_{m}J{\mathbf{C}}^{-1}\mathrm{Grad}m·{\mathbf{N}}_0\hfill \\ & =0\mathrm{on}\partial {\Omega }_0\times (0,t_{\mathrm{final}}].\hfill \end{array}$$

5. Numerical Tests

5.1. Discretisation and Implementation

We have employed mixed finite element methods for the spatial discretisation of systems (3a)–(5) and (28a)–(33d). We use quasi-uniform triangular/tetrahedral meshes in all cases. For the first system described in Section 2, since the linear poroelasticity equations are away from the nearly incompressibility limit (here the Poisson ratio is $\nu =0.4$), it suffices with using piecewise linear and continuous approximation of solid displacements (as well for all other unknowns). For the system proposed in Section 4, we employ the so-called Taylor–Hood mixed approximation of solid displacements and solid pressure composed by piecewise quadratic and continuous functions for each displacement component combined with piecewise linear and overall continuous functions for the solid pressure (this pair is known to satisfy an adequate inf-sup condition for the linearised hyperelasticity equations, see, e.g., [40]); and the fluid pressure as well as the chemical species concentrations are approximated with piecewise linear and continuous polynomials. For all problems, an outer fixed-point algorithm decouples the chemotaxis system from the poromechanics. An implicit-explicit method is employed to generate a full discretisation of the chemotaxis system and of the linear poroelasticity equations from Section 2, whereas for the poroelastic-growth system (28a)–(28c) with boundary conditions (32a)–(32b), an inner Newton–Raphson method is used to solve the corresponding nonlinear equations written in implicit form. All routines for the 2D tests have been implemented in the open-source finite element library FEniCS [41]. The solution of all nonlinear systems was carried out with Newton’s method using a tolerance of ${10}^{-7}$.on either the relative or absolute ${\ell }^2-$norm of the vector residuals, and the linear systems on each step were solved with the distributed direct solver MUMPS. For the 3D case we used the Firedrake library [42] (mainly due to its facility in handling block preconditioners). We provide details regarding the preconditioners in Section 5.3. The nonlinear tolerances are taken identical to the 2D case. For the linear solvers, we use a GMRES with the nested Schur preconditioner recently proposed in [8], with an absolute tolerance of ${10}^{-12}$ and a relative one of ${10}^{-2}$. The use of such a big relative tolerance yields a simplified inexact-Newton method [43] that results in more nonlinear iterations requiring much less time.

5.2. Mesh Independence Study

To experimentally assess the convergence of the proposed finite element method, we employ two tests. One study confirms that the Galerkin discretisation, when using the finite element spaces specified above, provides optimal (first-order) convergence with respect to mesh refinement (measured in the $H^1$-norm for all variables), where we simply use manufactured smooth solutions. The resulting convergence study is portrayed in

Table 1. Accuracy verification test using continuous and piecewise linear elements for the approximation of all variables. Error history (number of degrees of freedom, errors in the $H^1$-norm, and experimental convergence rates) with respect to smooth manufactured solutions computed at $t=10$. The symbol * indicates that the convergence rate is not computed for the coarsest level.

DoF	$\parallel \mathit{u}-{\mathit{u}}_{\mathit{h}}{\parallel }_1$	r	$\parallel \mathit{p}-{\mathit{p}}_{\mathit{h}}{\parallel }_1$	r	$\parallel \mathit{m}-{\mathit{m}}_{\mathit{h}}{\parallel }_1$	r	$\parallel \mathit{e}-{\mathit{e}}_{\mathit{h}}{\parallel }_1$	r	$\parallel \mathit{f}-{\mathit{f}}_{\mathit{h}}{\parallel }_1$	r	$\parallel \mathit{b}-{\mathit{b}}_{\mathit{h}}{\parallel }_1$	r
64	$1.0\times {10}^2$	*	$1.4\times {10}^3$	*	$1.0\times {10}^2$	*	$8.0\times {10}^0$	*	$2.8\times {10}^1$	*	$2.9\times {10}^1$	*
176	$3.0\times {10}^1$	1.77	$6.4\times {10}^2$	1.16	$2.8\times {10}^0$	1.88	$4.9\times {10}^0$	0.70	$1.9\times {10}^1$	0.67	$1.8\times {10}^1$	0.63
568	$1.1\times {10}^1$	1.19	$3.6\times {10}^2$	0.97	$1.1\times {10}^0$	1.36	$4.1\times {10}^0$	0.57	$1.4\times {10}^1$	0.53	$1.4\times {10}^1$	0.56
2024	$5.0\times {10}^0$	1.08	$1.7\times {10}^2$	1.19	$5.1\times {10}^{-1}$	1.09	$3.3\times {10}^0$	0.61	$8.5\times {10}^0$	0.73	$8.6\times {10}^0$	0.73
7624	$2.4\times {10}^0$	1.02	$8.3\times {10}^1$	1.11	$2.5\times {10}^{-1}$	1.02	$2.2\times {10}^0$	0.76	$4.4\times {10}^0$	0.96	$4.4\times {10}^0$	0.96
29,576	$1.3\times {10}^0$	0.98	$4.3\times {10}^1$	0.94	$1.3\times {10}^{-1}$	1.01	$1.3\times {10}^0$	0.82	$2.2\times {10}^0$	1.00	$2.2\times {10}^0$	0.99
116,488	$6.6\times {10}^{-1}$	0.98	$2.4\times {10}^1$	0.92	$6.3\times {10}^{-2}$	1.00	$6.5\times {10}^{-1}$	0.95	$1.1\times {10}^0$	1.00	$1.1\times {10}^0$	1.00

, where r denotes the experimental convergence rates computed as

$$r=log(\mathtt{err}/\widehat{\mathtt{err}}){[log(h/\widehat{h})]}^{-1},$$

where $\mathtt{err}$ and $\widehat{\mathtt{err}}$ denote errors produced on two consecutive meshes associated with mesh sizes h and $\widehat{h}$, respectively.

The second study determines the mesh resolution we employ to conduct the numerical simulations of growth in the linear and nonlinear regimes. For this we only run a test on the linear regime up to $t=100$ and simulating a single propagating wave. As an indicator we use a relative error (measured as the relative ${\ell }^2-$norm of the residual) with respect to a fine-mesh computation where the reference solution is computed on a mesh of approximately 300 K elements. The coarsest mesh I has 30 K elements and in the three cases we have used a graded mesh (locally refined near the edge $x=0$, where we expect the primation to initiate). The outcome of the mesh independence test is displayed in

Table 2. Mesh independence test. Approximate number of triangular elements and of degrees of freedom, relative computational time, and relative error with respect to a fine mesh computation.

Mesh	Cardinality	Total DoFs	CPU Time	Error
I	30 K	110 K	15%	9%
II	60 K	218 K	25%	5%
III	120 K	430 K	40%	4%

, from which we select mesh II for the subsequent 2D tests, since the relative error produced with this mesh is sufficiently small.

We also note that while the method does converge for much coarser meshes, grids as fine as meshes I, II, III are required to capture small features (for example, the spatial random distributions in the initial perturbations used, for example, in Section 5.4 below).

5.3. Efficient Preconditioners

Given that we are using a Newton solver, an efficient preconditioner is fundamental to obtain computational times that are feasible for the 3D case. The splitting strategy we are considering allows us to consider the poroelastic and chemotaxis problems separately, where each of them considers a different preconditioning strategy. In what follows we give the details of the preconditioners used for the solution of the tangent operator for each physics.

• Chemotaxis. This problem considers four similar building block physics, consisting essentially of three parabolic problems $(m,f,b)$ and one algebraic constraint $\left(e\right)$. We consider an additive block solver, meaning that we use a block-wise Jacobi preconditioner with only the diagonal block of each variable. At the block level, we consider the action of an AMG preconditioner for the parabolic problems and the action of a Jacobi preconditioner for the algebraic one.
• Poromechanics. The poromechanics block is more difficult, which is reflected in the complexity of the preconditioner under consideration. It is based on the block preconditioner proposed in [8] for large-strain poromechanics, where we use a lower Schur complement block factorisation with the fields $u$ and $(\psi ,p)$. For such a block we consider the action of an AMG preconditioner, whereas for the corresponding Schur complement block we consider instead a sparse representation given by an ad hoc extension of the fixed-stress splitting scheme proposed in [25]. For this, we add two stabilisation terms given by

  $${\int }_{{\Omega }_0}{\beta }_s\psi {\psi }^{*}+{\beta }_{f}p{p}^{*}dx,$$

  where we have denoted with ${(·)}^{*}$ the test function corresponding to each variable and have used the values ${\beta }_s={\beta }_f=0.1$. As this approach yields a sparse operator that approximates the Schur complement, we simply use the action of an AMG preconditioner based on it.

We highlight that the proposed preconditioner represents an improvement over [8], as we employ only preconditioner actions on each block, which allows us to use a GMRES linear solver instead of a flexible GMRES (which can be substantially more demanding in terms of computational cost). To assess the performance of the preconditioner, we consider only the first time instant of simulation, and report the nonlinear iterations, the maximum GMRES iterations per nonlinear iteration and the CPU times in

Table 3. Performance of the nonlinear solver used with the preconditioner detailed in Section 5.3. We show the total CPU time on the nonlinear solver on a single timestep, together with the maximum number of GMRES iterations in parenthesis.

DoFs	1 CPU	2 CPUs	4 CPUs	8 CPUs	12 CPUs	16 CPUs
3528	0.57 (10)	0.56 (10)	0.68 (10)	0.72 (10)	0.77 (10)	0.7 (10)
20,172	1.01 (10)	1.00 (9)	1.01 (10)	0.93 (10)	1.25 (10)	0.93 (10)
59,536	2.35 (11)	1.95 (9)	1.46 (10)	1.31 (9)	1.03 (10)	1.29 (10)
131,220	5.43 (11)	3.44 (10)	2.48 (10)	1.81 (10)	1.63 (10)	1.6 (10)
244,824	10.72 (11)	6.36 (11)	4.35 (11)	2.7 (10)	2.59 (10)	2.34 (10)
409,948	19.03 (12)	12.93 (12)	7.09 (12)	4.7 (12)	3.72 (12)	3.26 (12)
636,192	31.57 (12)	20.14 (12)	11.27 (12)	7.11 (12)	5.38 (12)	4.72 (12)
933,156	49.81 (12)	29.2 (12)	16.24 (12)	9.73 (12)	7.55 (12)	6.82 (12)
(a) Performance of chemotaxis solver.
DoFs	1 CPU	2 CPUs	4 CPUs	8 CPUs	12 CPUs	16 CPUs
16,893	2.5 (8)	1.8 (8)	1.8 (8)	2.05 (8)	1.91 (8)	2.17 (8)
108,501	12.13 (12)	7.27 (12)	4.8 (12)	3.11 (12)	2.8 (12)	2.89 (12)
337,229	39.71 (13)	22.62 (14)	12.31 (14)	7.57 (14)	6.6 (14)	5.72 (14)
765,477	98.6 (14)	52.06 (14)	28.45 (14)	16.5 (14)	12.57 (14)	10.79 (14)
1,455,645	202.15 (14)	102.76 (14)	56.51 (14)	32.32 (14)	23.83 (14)	19.58 (14)
2,470,133	350.59 (14)	192.52 (14)	99.84 (14)	55.04 (14)	40.24 (14)	33.75 (14)
3,871,341	598.98 (15)	324.14 (15)	168.48 (15)	89.02 (15)	63.27 (15)	54.29 (15)
5,721,669	869.06 (15)	462.3 (15)	242.9 (15)	130.09 (15)	98.15 (15)	78.02 (15)
(b) Performance of poromechanics solver.

.

5.4. Suppressed Solid Motion vs. Periodic Boundary Traction

In order to emphasise the effect of the stretch and deformation on the formation of patterns we compare the concentration of mesenchymal cells in two cases, one where the motion of the domain is suppressed (by setting $\tau =0$), and another scenario where the bottom of the domain, $\left(\mathsf{\Gamma }\right)$, is kept clamped, the top of the domain undergoes a (relatively small) periodic traction, and the vertical faces are kept with zero traction. In this case we employ the linear poroelastic model from Section 2, and a zero fluid pressure flux is imposed on $\mathsf{\Sigma }=\partial \Omega \backslash \mathsf{\Gamma }$. The spatial domain is a square $\Omega ={(0,20)}^2$ discretised into 58,482 triangular elements, and the model parameters that are modified with respect to those in Box 1 adopt the values

$$m_0=0.75,D_m=0.03,\tau ={\xi }_f=0.001.$$

The initial conditions correspond to uniformly distributed random perturbations of 10% with respect to the steady state $m_0,e_0,b_0$ and of 1% with respect to the steady state $f_0$. The traction is $t={(0,s_{0}sin(\pi t/\widehat{t}))}^{⊺}$, with magnitude $s_0=1500$ and period $\widehat{t}=320$. We employ a fixed time step $\Delta t=0.2$ and advance the system until $t=520$. The outcome of these tests is collected in Figure 5 and Figure 6, where we observed that slight modifications in the poromechanics result in quite drastic alterations of the clustering of mesenchymal cells. For instance, in Figure 6 we see that on the bottom edge of the domain (where displacements are zero during the whole simulation), the relatively large stretches applied elsewhere on the body do not break the formation of dotted-shaped patterns at early stages, but they are swept away as time advances.

5.5. Finite Growth in 2D

Next we maintain the domain and most model parameters as in the previous case and focus on the coupled model defined in Section 4. Determining the total first Piola–Kirchhoff stress from (27) and the growth specification (30) (with ${\delta }_2=0.15$, ${\delta }_3=0.045$), we modify also the boundary treatment and consider that ${\mathsf{\Gamma }}_0$ is the left edge of the domain whereas the remainder of $\partial {\Omega }_0$ is ${\mathsf{\Sigma }}_0$. On ${\mathsf{\Gamma }}_0$ we impose zero normal displacement $u·{\mathbf{N}}_0=0$ and a Robin boundary condition with constant ${\xi }_R=0.0001$ is set on ${\mathsf{\Sigma }}_0$. The shear modulus is $\mu =4$, the solid permeability is augmented to $\kappa =0.001$, and we take ${\xi }_f=0.3$. The fluid mass source is specified by (31) with ${\ell }_0=0.001$.

The panels in Figure 7 display transients of the growth of the initial domain and we plot the concentration of mesenchymal cells on the deformed configuration. We can observe that the patterns follow a different alignment than those formed in Figure 5, in particular disrupting the symmetry and the size distribution of the dotted-shaped patterns. The plots in the bottom row show the fluid and solid pressure together with the remaining chemical concentrations at the final time. Note that both solid and fluid pressures accumulate near ${\mathsf{\Sigma }}_0$. In addition, we note that the values of FGF at later times are negative since the last term on the left-hand side of (33c) represents a sink of f as the domain only grows (and $\mathrm{tr}\mathbf{C}-d$ is non-negative) throughout the simulation. We then change the sign of ${\xi }_f=-0.3$ and obtain positive FGF concentrations as well as similar patterns in all fields (but exhibiting a higher wave-number than in the previous case). The results from these tests are shown in Figure 8, which confirm a qualitative agreement with numerical tests from [6] and also with the actual patterns observed in nature.

5.6. Finite Growth in 3D

To close this section we proceed to extend the previous example to the 3D case. The undeformed domain is still a simple geometry ${\Omega }_0=(0,20)\times (0,20)\times (0,1)$, where to obtain an accurate solution we used 140 elements per side in the $X_1$ and $X_2$ directions, and instead considered only 8 in the $X_3$ direction. This results in roughly 10 million degrees of freedom. The sides at $X_1=0$ and $X_3=0$ constitute ${\mathsf{\Gamma }}_0$, where zero normal displacements are considered, and on the remainder of the boundary we set Robin conditions as above, but using ${\xi }_R=0.001$. The growth factors are considered as in (29) with the specification ${\delta }_1=1/52$, ${\delta }_2=0$, ${\delta }_3=1/104$. The results are shown in Figure 9, where we plot the evolution of the chemicals in the deformed configuration, and where the outline of the undeformed domain is also shown for visualisation purposes. We can see that as time progresses, the transversally isotropic growth depending only on time occurs mainly in the plane $X_1{X}_2$ as in the 2D case. Most interestingly, the mesenchymal concentration exhibits new complex structures. From the upper view, an arrow-like structure can be seen, where from higher concentrations appear more distant when going from right to left. From the slice view below, it can be seen how these distant concentrations are instead connected within the geometry. We stress that this type of patterns can only be described by the 3D model and currently there are no available results in the literature of mathematical and numerical methods to guide a qualitative/quantitative comparison. We also emphasise that the observed patterns are not necessarily representative of the actual primordial patterns—which are mainly composed by independent buds—and therefore this aspect still requires further investigation.

6. Concluding Remarks

We have proposed a model for describing the interaction between poroelastic deformations and chemotactical displacement of mesenchymal cells, FGF, BMP and epithelium. This new model is based on the theoretical framework developed in [6] for a rigid domain. In our first family of models, the coupling mechanisms include degradation of FGF due to local solid compression, advection of chemicals through solid velocity, and an active stress driven by mesenchymal cells. We have derived a set of conditions for Turing patterns to emerge, focusing on the two-way coupling and the main coupling variables. The analytical results on the stability of the steady solutions are confirmed by numerical solutions that concentrate on the regime of infinitesimal strains. Our numerical simulations reveal that the homogeneous oscillations predicted by the linear stability analysis are destabilised for sufficiently large values of the coupling parameters. We have also investigated an extended model that accounts for growth and hyperelastic deformations of the solid matrix, and the coupling mechanisms in this case consist in intrinsic growth functions depending on time and on the concentration of mesenchymal cells, as well as modification of the diffusion coefficients in the chemotaxis model due to changes from spatial to reference coordinates.

The linear stability analysis for the nonlinearly coupled system from Section 4 could be carried out extending our results from Section 3 following the ideas from [38]. As the main application is on skin patterning for feather development, we could also incorporate growth models specifically targeted for shells or thin plates as in [44], and concentrate on bi-layered structures extending the formulations in [45]. Viscous effects might be required in the rheological assumptions leading to the construction of stress, where relevant works in that context include [46,47]. Another interesting direction of development concerns performing further testing and optimisation of robust preconditioners that are required for the more computationally demanding 3D simulations. Finally, we mention that we are currently working on developing the well-posedness analysis of the fully nonlinear coupling of poroelasticity and chemotaxis models.