Numerical approximation of a mechanochemical interface model for skin appendage

We introduce a model for the mass transfer of molecular activators and inhibitors in two media separated by an interface, and study its interaction with the deformations exhibited by the two-layer skin tissue where they occur. The mathematical model results in a system of nonlinear advection-diffusion-reaction equations including cross-diffusion, and coupled with an interface elasticity problem. We propose a Galerkin method for the discretisation of the set of governing equations, involving also a suitable Newton linearisation, partitioned techniques, non-overlapping Schwarz alternating schemes, and high-order adaptive time stepping algorithms. The experimental accuracy and robustness of the proposed partitioned numerical methods is assessed, and some illustrating tests in 2D and 3D are provided to exemplify the coupling effects between the mechanical properties and the advection-diffusion-reaction interactions involving the two separate layers.


Introduction and modelling considerations
Reliable modelling of mechanochemical properties of composite materials is of key importance in a variety of engineering, material science and life science applications. These include tissue engineering, wound healing manipulation, bone fracture repair, design of photovoltaic devices, polymer adhesion, and many others. Of special interest to us are the two-way interactions between, on one hand, biological tissue deformation and, on the other hand, migration and proliferation of certain cell types during organismal development.
One application of such a general model is the study of the vertebrate skin appendage development and evolution, especially their diversity of forms (such as hairs, spines, feathers, and scales; [25,36,37,44]), and spatial organisations. Note that we have recently shown that all skin appendages are homologous structures developing from an anatomical placode and associated signalling molecules [17].
Turing proposed in his seminal work [59] that pattern formation in biological systems could be explained by a system of diffusion and interacting chemicals that he called 'morphogens'. Many subsequent studies [33,44,52,50] investigated biological patterning in the framework of such purely chemical reaction-diffusion models. However, as already suggested by Turing [59], mechanical aspects might also be relevant to biological patterning processes. During the last decade, biological experiments have clearly indicated that not only chemical (pattern of morphogen concentrations), but also mechanical (stresses and strains), parameters are key to skin development [1,3,25,42,58] and are synchronised during embryonic development [40].
Although the full detailed underlying mechanisms are still incompletely understood, authors have attempted to integrate both chemical and mechanical parameters in their mathematical models of skin patterning [41,43]. The various proposed mechanochemical models differ by the nature of the interactions between morphogen concentration and mechanics. For instance, while in [39,56] mechanical responses are directly triggered by chemical signalling, other models [38,39] propose that migration of cells directly generate local forces on the tissue. Such interactions between cells and morphogens can be complemented with chemical and mechanical exchanges between the two layers (dermis and epidermis) that constitutes the skin [15,51]. In one example [6], a three-chemical species model of limb bud development integrates growth velocity of the tissue normal to the limb surface and that directly depends on the local concentration of fibroblast growth factor. Other models consider the tissue growth as a liquid displacement [13], or aggregation of cells is the result of a chemical pre-pattern.
Here we focus on the computational modelling of mechanochemical interactions (inspired by [15] and [51]) based on simple reaction-diffusion systems governing the relations between, on one hand, the concentration of morphogenic proteins and cells densities and, on the other hand, mechanical stresses and strains. The mathematical description proposed in these two references is sufficiently general to be applied to multiple biological contexts (e.g., angiogenesis, fingerprint formation, cartilage condensation [41]), justifying the numerical algorithm presented here. Experimental evidence shows that the skin short-term and low-amplitude response to mechanical stimuli is that of an elastic solid [16]. We therefore assume that the two-layer domain of the skin (dermis and epidermis) can deform according to its inherent mechanical response, and that it is affected isotropically via an internal force that depends on the signalling factors. In addition, given the long spatio-temporal scale of the organismal developmental process, short-term fluctuations and inertial effects can probably be safely discarded to simplify the mathematical model [38]. We therefore adopt the two following fundamental assumptions: (i) the equilibrium of forces in the system is established by a quasi-static balance of linear momentum, and (ii) chemical species, governed by an advection-diffusion reaction (ADR) system, affect the medium deformation (through coupling functions). The chemical species' dynamics represents the spatio-temporal evolution of two morphogenic proteins and/or of cells that trigger and control the internal forces affecting the momentum balance. In turn, cell motion induces local and global traction forces on the elastic body, leading to substrate deformations [42]. We point out that cell motility depends on different factors that include kineses (increase in motile activity without any directional component), taxes (directed motion up or down a gradient), and guidances (motion directed by substratum cues, also called contact guidance) [39] that can also act locally or globally. Although long-range effects can be implemented through high-order spatial derivatives (appearing in the ADR or force conservation equations), we will here incorporate only short-range effects because our focus is on the numerical and algorithmic aspects of interface problems (such as the dermis-epidermis interactions in the skin).
As suggested elsewhere [15], here we distinguish between the two very different cellular dynamics of the dermis and the epidermis. Whereas we adopted a general logistic growth mitotic rate and non-negligible short range cell diffusion [40] in the dermis, we account for a zero production rate of cells and very small diffusivities in the equations governing epidermal variables [1,51]. This difference of treatment in the two tissues is justified by their respective structures: cells in epidermal tissue form sheets of adjoined interacting cells, while mesenchymal cells in the dermis are more mobile and separated by substantial extracellular matrix. However, as epidermal cells are not completely attached to each other and can rearrange topologically [1], we assume in both dermis and epidermis that all diagonal terms in the diffusivity matrices are strictly positive.
Regarding the interface problem between dermis and epidermis, we considered multiple methods for numerical approximation of linear elasticity using mixed finite elements: parallel and sequential Schwarz iterations using various interface conditions based on boundary element-finite couplings [23], hp-mortar [9], Nitsche [8], multiscale settings [11], or matching interface algorithms [62]. We finally chose a Schwarz iteration method with domain decomposition where the coupled set of mechanochemical equations is solved separately on each side of the interface, imposing mixed-type (Robin) boundary conditions, i.e., both fluxes and density values are equal in both tissues at their interface. A rich body of literature deals with choosing the weighting coefficients on the Robin interface conditions in order to optimise the convergence of the alternating Schwarz algorithms under different scenarios (see e.g. [12,21,26,55] and references therein). As the optimal weighting scheme tends to be problemdependent and might change substantially for diffusion, elasticity, or fluid-oriented problems, we will only describe general schemes and mention which specific weighting best suits the problems we investigate here.
Here, we define a family of finite element methods for the semi-discretisation of the problem, where the classical MINI-element [5] is used for the approximation of both displacement and pressure on each subdomain, whereas piecewise linear and continuous elements are used to discretise the vector of species concentration. Concerning the fully-discrete (spatial and time discretisation) scheme for the model problem, we apply an adaptive time-stepping strategy based on diagonally-implicit Runge-Kutta methods. This so-called TR-BDF2 method was introduced in the context of electronic circuit simulation [7] (see also [10,29,19] for further details) and has the advantage of using the same Newton iteration matrix to solve different temporal stages within a time step, a particularly interesting feature when addressing large systems of ODEs. The step size adaptivity exploits classical automatic control techniques used in the integration of so-called stiff systems [27,28], i.e., problems that tend to be numerically unstable.
This article is organised in the following manner. The remainder of this section contains preliminary notation and generalities to be used throughout the paper. Section 2 provides a detailed description of the differential equations constituting the model problem, with special emphasis on the coupling occurring at the interface. In Section 3, we introduce the Galerkin method, we present the adaptive time discretisation, and we specify algorithmic considerations. Section 4 contains a collection of computational tests focusing on the accuracy of the numerical schemes and on the suitability of the model equations to represent mechanochemical problems. We close with a summary and discuss potential extensions of the model in Section 5.
Throughout the paper, we let Ω Ă R d , d P t2, 3u denote a deformable body with polyhedral boundary BΩ, and denote by ν the outward unit normal vector on BΩ. Moreover, given a generic domain D P R d , d P t2, 3u, for two vectorial functions u, v P L 2 pDq d will denote their inner product as pu, vq D " ş D u¨v dx, whereas if S is a pd´1q´dimensional surface then their pairing will be written as xu, vy S " ş S u¨v ds. In addition, by P r pDq we will denote the space of polynomial functions of degree s ď r defined on the domain D, and I stands for the identity matrix in R dˆd . Finally, we will employ 0 to denote a generic null vector or the zero operator.

Set of governing equations
Linear elastic solids. We assume that the domain Ω is disjointly divided into two connected subdomains (with different material properties) denoted by Ω D and Ω E , respectively, whose intersection constitutes the interface between the regions and is denoted Σ " BΩ D X BΩ E . On BΩ E we define a part of the boundary related to the exposed surface, denoted as Γ (see Figure 2.1).
This heterogeneity of the skin (corresponding to the presence of an homogeneous dermis and an homogeneous epidermis, occupied by Ω D and Ω E , respectively) is accounted for by the spatiallydependent (but locally constant) Lamé coefficients Figure 2.1: Schematic representation of dermis and epidermis layers separated by the interface Σ and an exposed epidermic surface is denoted by Γ.
We also suppose that the two-layered elastic medium is nearly incompressible, so, in view of avoiding the well-known phenomenon of volumetric locking (inability of some discrete spaces to approximate the true displacement while satisfying incompressibility, see e.g. [9]), a mixed formulation is adopted in terms of displacements u˚and the pressure p˚"´λ˚div u˚, with˚P tD,Eu (although the interpretation of p as pressure holds exclusively in the incompressibility limit where λ˚{µ˚takes on quite large values). For all 0 ă t ă T , and given a body force f (whose specific form depends on the vector w; see equation (2.5) below) applied on each domain, we resolve the displacement / pressure pair pu˚ptq, p˚ptqq : Ω˚Ñ rR dˆR s on each subdomain, satisfying spring conditions at the surface boundary (Robin conditions), having clamped boundaries elsewhere, and assuming an adhesion condition at the interface (matching displacements and preservation of traction forces imply continuity of the medium across Σ). Then the global linear Navier-Lamé problem adopts the forḿ div`2µ˚εpu˚q´p˚I˘" f pw˚q in Ω˚ˆp0, T q, p˚`λ˚div u˚" 0 in Ω˚ˆp0, T q, where εpu˚q " 1 2 p∇u˚`∇u˚T q is the tensor of infinitesimal strains (symmetric gradient of displacement), σ˚" 2µ˚εpu˚q´p˚I is the Cauchy stress tensor, and α˚ą 0 is the Robin coefficient. The extension to non-homogeneous displacement conditions (u| BΩ ‰ 0) can be handled in a standard way by adding suitable liftings on the data. Below, we use the asterisk superscript˚P tD,Eu to denote quantities related to each subdomain, and we adopt the convention that the unit normal vector ν on the interface is pointing from Ω D to Ω E .
We stress that modelling the development of soft tissue and cells, in terms of mechanical dynamic behaviour is still an active area of research [2,61]. Biological tissues, depending on the nature of the constituent cells, present a large variability in their elastic properties (including viscous, nonlinear, or nonlocal effects) which can change substantially through spatio-temporal scales [35,45,61]. Some recent analyses were based on nonlinear St. Venant-Kirchhoff descriptions [35], including viscoelastic terms [2,39,57,59], and even fluid dynamics, as mentioned in [61]. A general consensus on how to represent mechanical properties of the skin tissue under morphogenesis is not yet available. Hence, our description will be restricted to classical linear elastic materials because it accurately represents small tissue deformations manifested during the first stage of skin appendage formation.
Advection-diffusion-reaction equations on an underlying deforming medium. Let us focus on the spatio-temporal interaction of the densities of m species w " pw 1 , . . . , w m q. We assume that all constituent species are equi-present at each spatial point (i.e., they are spatially continuous), so applying the Reynolds transport theorem to their mass conservation (assuming zero-flux boundary conditions at the surface Γ and continuity of fluxes and densities at the interface) leads to the system B t w˚`pB t u˚¨∇qw˚´divpM˚∇w˚q " Gpw˚q`gpu˚q in Ω˚ˆp0, T q, Here M˚P R mˆm is a (not necessarily diagonal, but positive definite) tensor of self-and crossdiffusion rates, G : Ω Ñ R m is a vectorial function containing the reaction kinetics of the system and representing the production and degradation of species concentration on each subdomain Ω D , Ω E ; and the term g encodes an additional coupling contribution (detailed in equation (2.5) below). The system is endowed with adequate initial conditions for w˚, and the general form (2.2) accommodates either dimensional or dimensionless descriptions.
Here,we focus on two levels of model complexity. Firstly, we define a simple two-variable activatorinhibitor model with quadratic reaction terms due to Gierer-Meinhardt (cf. [24]): where w1 is the short-range autocatalytic substance and w2 is the long-range antagonist; and the ρ i 's are positive model parameters.
Second, we introduce a more general model of skin appendages development were w1 is the cell density, w2 is the density of the extracellular matrix, and w3 , w4 are the concentrations of two morphogen signalling chemical compounds. All quantities are assumed to diffuse, so none of the diagonal terms in M˚is zero (nevertheless, self-diffusion for cell concentration in the epidermis and the matrix density in the two tissues, will be very small in comparison to the other values in the diagonal of M E ), and we discard higher order terms (derivative with order higher than two) representing haptotaxis and non-local effects in the cell migration fluxes. The reactions consist in a logistic relation for the mitotic rate, a zero rate of matrix material secretion, and the production, as well as degradation, of morphogens: where r 0 is the maximum cell density, and r i 's are positive model constants. The specifications in (2.4) have been adapted from the skin pattern formation model proposed in [15]. Similar models can be found in the context of cell migration, angiogenesis, follicle differentiation, or mesenchymal morphogenesis (see e.g. [1,34,37,43,60]).
In the specific context of skin appendage formation, signalling pathways control cell movements that eventually lead to mechanical tissue deformation due to cell aggregation [1,25,53]. These interactions can be represented by diverse local or non-local processes including advection, cytoskeleton regulation, chemical gradients, active stress, strain dependent diffusivity, and many others [40,38]. In the present model we opt for including only local variations of the ADR and mechanical variables. This simplification rules out some phenomena such as inter-cellular force transmission through cytoskeleton, nonlinear interactions (e.g., cells sensing, through filopodia, long-range variations), or anisotropy effects, which can be of special importance in some biological contexts [15,34]. Nonetheless, as suggested in [38], local models such as the one we use here, may suffice to reproduce accurately a variety of mechanical processes without the need for additional constitutive relations.
The PDE-based model (2.1)-(2.2) (written in Eulerian form and solving for displacements, solid pressure, and molecular variables) assumes that changes in the chemical concentrations do not affect the mechanical properties of the solid, nor the Lamé constants; the species only act as forcing terms on the momentum balance, i.e., the intensity of the deformation will depend locally on the gradient of the species concentrations. Conversely, we here suppose that the body deformation affects the species dynamics by means of advection and through the term gpu˚q on the right-hand side in (2.2), carrying local information about the medium dilation. Therefore, we set where cf , cg are model parameters. Note that the form of f is equivalent to considering a stress component depending on the concentrations of some wi , that is σ˚" σe ff`wi I, where σe ff is an effective elastic stress. Similar coupling conditions can be found in [22,31,49].
The unique solvability of (2.6) has not yet been addressed in the literature, however related results concerning the coupling of linearised mechanics and reaction-diffusion systems have been recently addressed in the context of cardiac and plant cell biomechanics [4,46].
Treating the interface: a non-overlapping continuous Schwarz method. As mentioned previously, the treatment of the interface (here, between the epidermis and dermis) is an important issue for interacting media problems and this paragraph will focus on its realisation. Regarding the geometric separation between the problems defined on each side of the interface, we remark that (2.6) is valid in particular if we consider a solution globally defined in Ω. Nevertheless, in the discrete case we aim at solving two separate problems on each subdomain and therefore we require to impose jump conditions in an adequate manner. We adopt a non-overlapping Schwarz iteration based on Robin-Robin type interface conditions for the ADR system where Rp¨q denotes the Robin transmission operator, and K D , K E are nonnegative acceleration constants satisfying K D`K E ą 0 (see, for instance [32] and [18,47]). On the other hand, the adhesion interface condition featured by the elasticity problem is incorporated through the following relation where Sp¨,¨q is a transmission operator associated to the elasticity problem, and J˚ą 0 is another acceleration constant. The classical iterative Schwarz algorithm (still continuous and still nonlinear) then reads: For t ą 0 and for " 0, . . ., find pu˚, `1 , p˚, `1 , w˚, `1 q P V˚ˆQ˚ˆZ˚such that J˚on Σ, 0 otherwise, and the superscript : P tD,Eu denotes quantities in the subdomain opposite to˚. Here we have dropped the explicit time-dependence whenever clear from the context.
Note that the order in which these blocks are solved is arbitrary. We also note that the blocks (3.2),(3.5) still constitute nonlinear problems, that we solve through inner Newton iterations. Alternatively, one can commute the domain decomposition method and the nonlinear solvers, in such a way that the nonlinearly coupled set of equations can be solved monolithically on each subdomain. Further details on the choice we adopt here, are provided in Algorithm 1, below. Finally we remark that the sequential steps (3.2)-(3.5) can be regarded as a block Gauss-Seidel iteration applied to the fully monolithic differential-algebraic system.
For the discretisation of the elasticity equations in mixed form, we require the discrete spaces in (3.1) to be inf-sup stable. In our implementation we choose the so-called MINI elements [5], characterised by Vh :" tv P V˚: v| K P rP 1 pKq`BpKqs d , K P Th u, Qh :" tq P Q˚: where BpKq is the space of cubic bubble functions defined locally on an element K; and we employ the space of conforming Lagrangian finite elements for the approximation of the species concentrations Zh :" tz P Z˚: z| K P P 1 pKq m , K P Th u.
The discretisation (3.2)-(3.5) leads to the following system of equations M˚9 w˚, l`1`Ců w˚, l`1`D w˚, l`1´G˚p w˚, l`1 q " 0, where the superscript ; may be replaced by the appropriate index associated to the Schwarz algorithm. We readily observe that only (3.6) is a nonlinear ODE system, whereas the matrix form of the elasticity sub-problem consists of a linear system. The subscript u in Ců emphasises that the matrix is directly dependent on displacements.

Adaptive timestep scheme
The family of fully-discrete methods that we use here is based on adaptive, embedded Runge-Kutta schemes of order 2 and 3 (RK23). For a generic ODE system d dt v h " L h pv h , ϕ h ptqq, we recall that the main step in the classical RK23 scheme reads with s " 3 stages, and where the two approximations v h , p v h are of order Oppq and Opp pq, respectively. As in classical embedded RK methods, here the coefficients ta ij , c j u coincide in both approximations; however the coefficients b j differ in order to attain the expected convergence orders (in general one chooses p p " p˘1, cf. [27]). Even when the stability properties of the two schemes do not match, embedded RK methods provide a rather non-expensive estimator for the local error (of the less accurate solution) which is then directly used in the control of the timestep [27,28]. The specific description of the controller used here is presented in Algorithm 2, below. The RK23 scheme is used together with the TR-BDF2 method (see [7]), which is a particular class of Diagonally Implicit Runge-Kutta (DIRK) methods (that is, a ij " 0 for i ă j, and at least one a ii ‰ 0) [27]. The scheme is characterised by the following Butcher's tableau reported in Table 3 and γ " 2´?2 2 . As established in [29], the A-stability, the stiffly accuracy (a sj " b j , for j " 1, . . . , s), as well as the strong S-stability of the method are guaranteed by the coefficients used in Table 3.1. An appealing advantage of this scheme is that, for each step, only two implicit stages are required (where one can reuse the same Newton iteration matrix) while one stage is explicit. In this regard, the method is closely related to singly DIRK (SDIRK) methods (where all diagonal a ii coefficients are equal and the first stage is explicit) [29]. In addition, TR-BDF2 falls in the class of First-Same-As-Last (FSAL) RK schemes, i.e. methods where the current explicit stage coincides with the last stage of the previous step, therefore decreasing the overall computational cost.
These problems should be solved through Newton iterations before moving to the next timestep. Alternatively, and as suggested in [28], one can apply the change of variables v . We follow this approach and modify the superscripts p1q Ñ n, p2q Ñ n`γ and p3q Ñ new, to obtain the following full discretisation of the ODE system: w˚, n`γ " γw˚, n`γ z`γ w˚, n z`w˚, n , w˚, new " γw˚, new z`b2 w˚, n`γ z`p 1´b 2´γ qw˚, n z`w˚, n , with w˚, n`γ z " pM˚q´1∆t n "´Ců w˚, n`γ´D w˚, n`γ`G˚p w˚, n`γ q ‰ , w˚, new z " pM˚q´1∆t n r´Cůw˚, new´D w˚, new`G˚p w˚, new qs .
This nonlinear system is solved by means of Newton steps as follows wz ,k`1 " wz ,k`δk`1 wz , We finally point out that an additional reduction of computational cost occurs due to the fact that the local error err " p

Algorithmic details
Overall solution algorithm. The global numerical scheme to solve the mechanochemical model is summarised in Alg. 1. After assigning suitable initial conditions u˚, 0 " 0, p˚, 0 " 0, w˚, 0 " w0 , projected to the corresponding finite element spaces, we use the non-overlapping Schwarz method (3.2)-(3.5) to produce a solution until time T f . Two important points stand out from this master algorithm. First, only the ADR equations are solved through the embedded RK method within the TR-BDF2 method (see Alg. 1, Line 8, detailed in Alg. 2), whereas the elasticity subproblem is solved by an implicit method (cf. Line 11 in Alg. 1). Furthermore, the velocity in the advection term is approximated from the solid displacement using a backward Euler scheme (cf. Alg. 1, Line 10 and 12). Second, the computation performed in Alg. 1, Line 5 and 18, represents an alternative way of performing the explicit first stage of the RK scheme [29]. Such a rescaling must be executed at each modification of the timestep, so that one starts with an adequate value for w z when entering Line 8 of Alg. 1. An additional benefit of the rescaling is that it improves the stability of the overall stiff problem, and at the same time avoids unnecessary function evaluations.
Adaptive Runge-Kutta method. The RK procedure is detailed in Alg. 2 and it features a similar structure as the ODE solver ode23tb in MATLAB. Again, we highlight two major components: the one-step RK resolution and the adaptive setting of the timesteps. Line 2 in Alg. 2 represents stage 2 of the RK method, i.e., the first implicit stage (S1). As the scheme is implicit, we initialise the tangent system associated to the Newton iteration. These matrices will actually be assembled only if the Boolean mjcontrol is true, and this switch is described in Alg. 4, below. Should mjcontrol be false, the mass-Jacobian blocks are evaluated at each Newton step (see Alg. 2, specifically Line 8). We also stress that the initial guess for the rescaled derivative w n`γ z can be simply taken as w n z , since extrapolating the derivative of the interpolant from the previous step does not improve accuracy nor efficiency. Next, when updating the solutions we also change the Robin transmission conditions so that the Schwarz algorithm is performed at the same time (see Line 14 in Alg. 2). In order to retrieve w new z,k`1 and w new k`1 , a second Newton step is necessary within the implicit stage (S2), the particularity of which is the evaluation of the initial guess w new z,0 , by extrapolation of the derivatives of the cubic Hermite interpolant from t n to t n`γ [29].
The algorithm also focuses on the automatic determination of each timestep, in view of the stability of the time integration of stiff problems. The process is designed to reject the overall step if one of the two implicit stages fails or if the local error is too large. The potential rejection at an implicit stage is monitored by the Boolean variables rejectS1 and rejectS2, which are switched on according to the value of the residual in the Newton iteration described in Alg. 3. How the value is actually computed may affect the performance of the algorithm (see e.g. [54]). For a generic residual vector v of length M , here we adopt the rescaled maximum norm v~" max i"1,...,M pv i¨s c´1 i q, Algorithm 1 Overall solution algorithm for the mechanochemical problem 1: Initialise the solution vectors u˚, 0 , p˚, 0 , w˚, 0 and define w fi ∆t min " |t|, ∆t old " ∆t, ∆t " minp∆t max , maxp∆t min , ∆t old qq 4: if ∆t ‰ ∆t old then set n " n`1, tn`1 " tnew, w n`1 " w new , w n`1 z " w new z , and |w n`1 |2 " |w new |2

19:
needNewMJ Ð true 20: end if 21: end while where sc i :" maxpsc i , |w k`1,i |q is the local maximum of the solution vector at the k-th Newton step, and it is reinitialised as sc i " maxp|w n | 2 , ηq, where η " ATOL RTOL . This process is summarised in Line 1 from Alg. 2.
We recall (see Line 1 of Alg. 3) that the residual in the Newton step is simply the norm of the current increment δ k`1 wz . The loop can break under different cases. For instance, if the residual is smaller than a given threshold, then Line 2 terminates the step. Otherwise, the linear convergence of the overall scheme (see e.g. [28]) implies that one can evaluate the iteration error using the residual together with the inequality }δ k`1 w z } ď θ}δ k w z }, with θ ă 1. This experimental rate is reset to zero at each mass-stiffness assembly so that it depends on the Newton iteration matrix [48]. The rate is then estimated using θ 1´θ res (see Line 20) with two successive iterations and the current error. We consider that the algorithm has converged if the error is smaller than κ N times a fixed tolerance. A given stage is therefore rejected if at least one of the following situations occurs: (a) the convergence rate is larger than 1 (cf. Line 16); (b) the maximum number of iterations K " 10 is reached (cf. [28] and Line 23); and (c) For any k, we have κ N TOL N ă epδ k`1 w z , p scqθ K´k (cf. Line 26). Point (c) represents a rough estimate of the iteration error expected after K´1 iterations.
Should a stage be rejected, one moves to Line 25 in Alg. 2 and reduces the timestep ∆t by a factor f ac S1S2 (independently of the mass matrix assembly). Then, Line 31 checks whether both mass and Jacobian matrices have been updated. Finally, Line 36 verifies that the local error estimation remains below a prescribed tolerance. Based on the estimators for the scaled derivatives we proceed to compute the local asymptotic error e 1 w as a result from the linear system stated in Line 42 of Alg. 2. This step is necessary as the boundedness of the raw estimator is not always guaranteed [28].
Next, we select the maximum error between these two estimators as an approximation of the timestep error. In contrast with [28], we do not observe here the so-called hump phenomenon intrinsic to automatic step methods, so we do not require an additional correction step. We only update the timestep in an automatic fashion (thanks to the sequence of rejection procedures), capturing correctly the dynamics of the integrated system. When all errors are lower than the corresponding tolerance, the timestep is accepted and the interface elasticity problem is solved in Line 11 from Alg. 1. We stress that the algorithm will stop if ∆t ď ∆t min . Algorithm 2 Adaptive Runge-Kutta time stepping: procedure TR-BDF2(w n ,w n z ,∆t,β ∆t ) 1: compute sci " maxp|w n |2, ηq 2: assemble M and J using MJcontrol(w n ,w n z ,∆t,β ∆t ) Ź See Algorithm 4 3: set tn`γ " t`2γ∆t, w n`γ 0 " w n`2 γw n z , w n`γ z,0 " w n z Ź Stage 2 4: rejectS1 Ð false 5: set resS 2 ,min " 100 ~w n`γ 06 : for k " 0, . . . , K do solve δ n`γ k`1 wz " p 1 ∆t M˚´γJ˚q´1g˚, n`γ k 13: compute w˚, n`γ z,k`1 " w˚, n`γ z,k`δ n`γ k`1 wz , w˚, n`γ k`1 " w˚, n`γ k`γ δ n`γ k`1 wz compute sci " maxp|w n`γ |2, sciq 21: w new z,0 " p31w n z`p32 w n`γ z`p33 w n`γ´p if mjcontrol then 27: if Jcurrent^Mcurrent then 28: set ∆t old " ∆t, ∆t " maxpf acS1S2∆t old , ∆t min q, w n z " ∆t ∆t old w n z 29: needNewMJ Ð true if err ą RTOL then 44: ∆t old " ∆t, ∆t " maxˆ∆t old maxˆf acmin, f ac´R TOL err¯k I q, ∆t min˙˙, w n z " ∆t ∆t old w n z 45: needNewMJ Ð true 46: end if 47: end if Alg. 4 permits to control the mass-stiffness assembly. This procedure consists of a number of controllers, including the aforementioned mjcontrol (see a list in Table 3.2). The state modification of these variables depends on step-size revaluation, which itself depends on the rejected stage of the RK method. Notice that Mcurrent and Jcurrent control whether the mass and Jacobian matrices are evaluated at each new timestep (if this timestep is rejected at least once). For models such as (2.2), the label Mcurrent is always true. However if we want to solve the interface problem using e.g. a level-set approach, then the mass matrix needs to be recomputed and the controller can be Algorithm 3 Newton stopping criterion: procedure NewtonConv(w k`1 ,δ k`1 w z ,res min ,β ∆t ) 1: compute p sci " maxpsci, |w k`1,i |q and res k`1 p p sciq "~δ k`1 w z2 : if res k`1 ď resmin then

Numerical tests
This section contains a set of examples assessing the experimental spatio-temporal convergence of our primal-mixed method, and also illustrates its use in two mechanochemical interface problems. The implementation is carried out with the open-source finite element library FreeFem++ [20], and all sparse linear systems are solved with distributed direct methods (MUMPS and UMFPACK).  Table 3.2: Set of Booleans controlling the evaluation and assembly of mass and stiffness matrices. Example 1: Experimental accuracy against manufactured solutions. Let us consider a rectangular domain Ω " p0, 1qˆp0, 1.4q divided into the dermis Ω D " p0, 1q 2 and epidermis Ω E " ΩzΩ D subdomains. The initial position of the interface is characterised by the segment x P p0, 1q, y " 1, and the external surface is x P p0, 1q, y " 1.4. The main goal of our first series of tests is to analyse the convergence of the spatial discretisation. We therefore focus on a stationary problem where no advection occurs in the morphogen model. We construct smooth exact solutions, defined as follows: w "ˆ1´c osp2πxq sinp3πyq 1`1 2 cosp2πxq sinp3πyq˙,ũ "ˆx p1´xq cospπxq sinp2πyq sinpπxq cospπyqy 2 p1´yq˙, (4.1) and set u D "ũ| Ω D , u E "ũ| Ω E , w D "w| Ω D , w E "w| Ω E , and p˚"´λ˚divũ˚. The functions in (4.1) are used to construct the Dirichlet datum for the displacements on BΩzΓ as well as the flux conditions imposed on the morphogens at the boundary BΩ. Suitable non-homogeneous forcing terms are also set using these closed-form solutions. We choose the two-species Gierer-Meinhardt reaction model (2.3), and adopt the coupling terms as in (2.5). Model parameters are chosen as in Table 4.1.
We proceed to generate a sequence of successively refined unstructured triangular meshes for both subdomains, and compute errors between exact and approximate numerical solutions (recall that we are solving for the steady version of the coupled problem) on each refinement level. For a generic scalar or vector field s˚we will denote errors and convergence rates as follows e 0 ps˚q :" }s˚´sh} 0,Ω˚, e 1 ps˚q :" }s˚´sh} 1,Ω˚, r k ps˚q " logpe k ps˚q{ p e k ps˚qq logph{ p hq , k " 0, 1, where e k and p e k denote errors computed on two consecutive meshes of sizes h and p h, respectively. The outcome of the convergence test (using the MINI-element for the approximation of displacement and pressure, and Lagrangian finite elements for the species concentration) is presented in the first part of Table 4.2, which indicates that the method converges optimally. Four Newton iterations are sufficient for each mesh to reach the tolerance of TOL N " 1e-5, and the stopping criterion is based on the 8´n orm of the total residual. The time convergence is assessed by taking a fine mesh of size h " 0.0013 and computing approximate solutions to the transient coupled problem (2.1),(2.2) compared against the exact displacement and pressure from (4.1) whereas the species concentrations are now w "w expp´ktq,  with k " logp1{2q. Errors and convergence rates are measured and denoted as follows e ∆t ps˚q :" N ÿ n"1 }s˚pt n q´shpt n q} 0,Ω˚, r ∆t ps˚q " logpe ∆t ps˚q{y e ∆t ps˚qq logp∆t{ x ∆tq , where e ∆t , y e ∆t denote errors computed using two consecutive timesteps ∆t and x ∆t, respectively. The second block in Table 4.2 shows an asymptotic second order convergence, as expected from the use of the TR-BDF2 time stepping algorithm (see also [10,19]  Example 2: Cross-diffusion effects in a two-dimensional morphogen interface model. In this example we focus our attention on the effect of cross-diffusion Mi j ‰ 0 @i, j in (2.2). We consider again the Gierer-Meinhardt kinetics (2.3), now decoupled from the linear elasticity equations (2.1). We impose the same reaction parameters on both sides of the interface ρ 0 " 0, ρ 1 " ρ 2 " ρ 4 " ρ 5 " 1, ρ 3 " 0.35, we set the diffusion coefficients M 11 " M 1 " 1, M 22 " M 2 " 30, and choose the acceleration transmission constants K D " K E " 1. The domains of interest are the blocks Ω D " p0, 50q 2 and Ω E " p0, 50qˆp50, 75q, which we discretise into unstructured meshes with 16360 and 8032 triangular elements, respectively. Initial concentrations are set according tõ where ηpxq is a uniformly distributed random field with variance 10´3. We employ the adaptive time-stepping method described in Algorithm 2 with controller parameters defined as in Table 4.3.
We simulate the process until T " 2000, and display the resulting numerical solutions on  In Figure 4.2 we observe that the adaptive time stepping behaves differently depending on whether one reuses or not the Newton iteration matrices over the implicit stages (this feature is turned on/off by the parameter mjcontrol). For each case, we observe an increase of ∆t as the system reaches a steady state. The spatial patterns generated by the two approaches practically coincide and therefore are not shown. However, we observe that the case mjcontrol " true does not show a marked increase of the timestep when the system reaches stationary patterns. Additionally, we observe that reusing the iteration matrix leads to a slight increase in the number of Newton iterations in the S2 stage.
In addition, we test our algorithm for both linear and nonlinear cross-diffusion terms defined as where Mi j , η˚, w k ij (i, j, k " t1, 2u) are real constants. For these simulations we impose M1 1 " M1 , M2 2 " M2 . As illustrated in Figure 4.3, the final patterns adopt different structures according to the cross-diffusion matrix employed. While the specific nonlinear cross-diffusion matrix modifies the size of w i -concentration spots without changing its structure, the linear one (and depending on the parameter values) can change from spots to striped patterns. Furthermore, nonlinear cross-diffusion induces a more visible impact on the distance between high concentration regions, when compared to the linear case. Space-parameter plots illustrate this phenomenon in Figure 4.4. For both crossdiffusion models, we vary spatially the values rM1 2 , M2 1 s P r0, 1sˆr0, 15s, and η˚, w1  " η˚, w2

22
" η 2 from 0 to 1; always checking that the cross-diffusion matrices remain positive definite. The linear cross-diffusion particularly affects the type of patterns that are generated but does not significantly change length scale of the pattern. Nonlinear diffusion, on the other hand, affects directly the pattern length scale without altering the general pattern structure produced by the system. These observations are confirmed in 3D domains as well (cf. 4.5 ).
Example 3: A simplified activator-inhibitor model on an elastic 2D domain. With the objective of testing the time-adaptive algorithms now applied to the mechanochemical coupled problem, we add the elasticity components to the reaction kinetics (2.3) of Example 2 (with the same parameters and mesh). The Young moduli and Poisson ratios defining the material properties of the medium are now E D " 1000, E E " 250, ν D " 0.475, ν E " 0.3. Displacements are prescribed on the boundary BΩzΓ according to u " p 1 2 cosp5πpx´yq{2q, 3 4 sinp5πpx`yq{2qq T . Additionally, we impose c D f " 150, c E f " 20 and cg " 1. The medium is initially at rest and stress-free (zero displacement and pressure on the whole domain) and we re-use the initial conditions of the morphogen concentrations defined in the previous test. We use Algorithm 2 with controller parameters defined again as in Table  4.3, except for A TOL " 10´6 and TOL N " R TOL " 10´4. The process is simulated until T " 2000.  Figure 4.1 (in terms of magnitude and spatial structure). For the layer with smaller Poisson ratio, the patterns tend to homogenise and eventually disappear. This phenomenon highlights the effects of mechanical properties on the morphogen dynamics, and we also note that deformations are larger near the interface and close to the bottom boundary of Ω D .
are responsible for considerable variations on the species concentration spatial distribution as well as on the deformation profiles produced by the model. For instance, we observe that the size and shape of morphogens concentration patterns induce differences on the deformation of the interface Σ. Conversely, differences on the mechanical response of each subdomain also directly impact on the chemical species distribution, at the interface. Surprisingly, varying the coupling parameters does not substantially impact the length scale but rather make these patterns to disappear altogether.
We also inspected the effect of the acceleration constants of the ADR interface conditions together with the degree of material heterogeneity (how different the Lamé parameters are in each subdomain) on the timestep evolution and on the Newton iteration count. For these simulations we fix mjcontrol = true in order to inspect also the efficiency on the update of the Jacobian matrix throughout the simulation. Some of these results are collected in Figure 4.8. The value N SX N ewt represents the number of iterations corresponding to accepted timesteps. For this case we observe that ∆t is of the order of 10´3 and progressively increases with the simulation time and the consolidation of the spatial patterns. The algorithm readily detects when the process approaches a steady state, leading to a fast increase of the timestep. However, as show in the case using ν˚" 0.3, this behaviour depends on the choice of model parameters and oscillations might appear, indicating a suboptimal adaptive time stepping. Moreover, changing the value of the acceleration constants may also affect the simulation time, the deformed interface, and the pattern structures. In all our experiments, we observe that accepted steps Example 4: Simulation of the full-3D appendage model. We finally turn to the simulation of the mechanochemical interaction in the context of skin appendage using a 3D slab as computational domain, and focusing on the reaction model proposed in (2.4). Model parameters are chosen according to Table 4.4. The domains of interest are the blocks Ω D " p0, 50q 3 and Ω E " p0, 50q 2ˆp 50, 75q, which we discretise into unstructured meshes with 104360 and 60067 tetrahedral elements, respectively. Initial concentrations are set according to steady state solutions obtain from (2.4), perturbed with a uniformly distributed random variable with variance 10´3. The adaptive method described in Algorithm 2 is employed with controller parameters as in Table 4.3, except for A TOL " 10´6 and TOL N " R TOL " 10´3. The process is run until T " 300.
The results are presented in Figure 4.9, showing a clear shift in species (cell, matrix and morphogens) concentration patterns across the interface, with the exception of species w2 (explained by the similarity in the reaction description (2.4)). The uniform Lamé constants on both domains and the Robin conditions imply that deformation patterns are of comparable magnitude throughout Ω, and the zones of high deformations are concentrated near the domain boundary. Different mechanochemical patterns can be produced by changing the coupling and reaction parameters, and a specific study on linear and nonlinear stability will be carried out in a forthcoming contribution. At this stage we can already observe that for w1 , for example, only the dermis species possesses a production term (with logistic growth) while the dynamics on the epidermis is described by a convection-diffusion process. This explains why the concentration variation in the dermis is much smaller than in the epidermis layer. For w3 and w4 we observe a similar behaviour, but having a larger variation on the z-axis for a specific layer (Ω D for w 4 , Ω E for w 3 ). This is due to some chemical species being produced only on a specific layer (see (2.4)) and the relevant concentration decreases rapidly, due to metabolism . Linear cross-diffusion (top row) using pM1 2 , M2 1 q " p0.5, 0q, p0, 0.5q and p0.35, 0.05q (left, centre, right, respectively). The bottom panels show the case of nonlinear cross-diffusion using pη˚, w1 11 , η˚, w2 11 , η˚, w1 22 , η˚, w2 22 q " p10´1, 10´1, 10´1, 10´1q, p1, 0, 0, 1q and p1, 0, 1, 0q (left, centre, right, respectively). degradation, as one moves away from that layer.

Concluding remarks
We have introduced a mechanochemical model for the simulation of basic processes related to skin appendage patterning mechanisms. The two-way coupling between linear elasticity and advectiondiffusion-reaction systems is achieved through a mechanical term depending on the gradient of the chemical species, and a pressure-dependent source term representing production of species. A weak formulation in a multidomain setting was introduced, and we presented a partitioned fixed-point and Schwarz algorithm to decouple the system into four principal blocks (the two domains, dermis and epidermis, being treated for advection-diffusion-reaction and mechanics). The space discretisation uses MINI-elements for the approximation of the displacement-pressure pair, and piecewise linear and continuous Lagrange elements for the species concentrations. A trapezoidal BDF2 method combined with adaptive Runge-Kutta schemes is used as a time advancing strategy, and algorithmic details are provided. The convergence properties of the proposed methods were studied in detail, and a set of numerical tests addresses pattern generation and the deformation of the interface.
Anomalous diffusivity (e.g., fractional, stress-or strain-dependent diffusion [14]) was not considered here. Further work is also needed on the detailed linear stability analysis associated to the coupling mechanisms for the specific models at hand, as well as on the application of the developed methods in the simulation of more complex phenomena. Future work will also address extension of our model to non-linear elasticity and viscoelasticity and theoretical considerations on the convergence properties of partitioned and monolithic schemes.