# Numerical approximation of a mechanochemical
interface model for skin appendage††thanks: This work was partially supported
by grants from the Swiss National Science Foundation (FNSNF, grants
31003A_140785 and SINERGIA CRSII3_132430), by the SystemsX.ch initiative
(project Epi-PhysX); and by the Engineering and Physical Sciences Research
Council (EPSRC, research grant EP/R00207X/1).

Luis Miguel De Oliveira Vilaca, Michel C. Milinkovitch†, Ricardo Ruiz-Baier
Laboratory of Artificial & Natural Evolution (LANE), Department of Genetics
and Evolution, University of Geneva, 4 Boulevard d’Yvoy, 1205 Geneva,
Switzerland; and SIB Swiss Institute of Bioinformatics, Geneva, Switzerland.
Email:
{LuisMiguel.DeOliveiraVilaca,Michel.Milinkovitch<EMAIL_ADDRESS>author. Mathematical Institute, Oxford University, Andrew Wiles Building,
Woodstock Road, OX2 6GG Oxford, UK. Email<EMAIL_ADDRESS>
###### Abstract

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.

Key words: Elasticity-diffusion problem; Interface coupling; Pattern
formation; Skin appendage modelling; Finite element methods; Adaptive time
stepping.

Mathematics subject classifications (2000): 92C15, 35Q92, 65M60, 74S05, 35K57,
74A50.

##  1 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 problem-dependent 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 $\Omega\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$
denote a deformable body with polyhedral boundary $\partial\Omega$, and denote
by $\boldsymbol{\nu}$ the outward unit normal vector on $\partial\Omega$.
Moreover, given a generic domain $\mathcal{D}\in\mathbb{R}^{d}$,
$d\in\\{2,3\\}$, for two vectorial functions $\boldsymbol{u},\boldsymbol{v}\in
L^{2}(\mathcal{D})^{d}$ will denote their inner product as
$(\boldsymbol{u},\boldsymbol{v})_{\mathcal{D}}=\int_{\mathcal{D}}\boldsymbol{u}\cdot\boldsymbol{v}\,\mbox{d}\boldsymbol{x}$,
whereas if $\mathcal{S}$ is a $(d-1)-$dimensional surface then their pairing
will be written as
$\langle\boldsymbol{u},\boldsymbol{v}\rangle_{\mathcal{S}}=\int_{\mathcal{S}}\boldsymbol{u}\cdot\boldsymbol{v}\,\mbox{d}s$.
In addition, by $\mathbb{P}_{r}(\mathcal{D})$ we will denote the space of
polynomial functions of degree $s\leq r$ defined on the domain $\mathcal{D}$,
and $\mathbf{I}$ stands for the identity matrix in $\mathbb{R}^{d\times d}$.
Finally, we will employ $\boldsymbol{0}$ to denote a generic null vector or
the zero operator.

##  2 Set of governing equations

#### Linear elastic solids.

We assume that the domain $\Omega$ is disjointly divided into two connected
subdomains (with different material properties) denoted by
$\Omega^{\mathrm{D}}$ and $\Omega^{\mathrm{E}}$, respectively, whose
intersection constitutes the interface between the regions and is denoted
$\Sigma=\partial\Omega^{\mathrm{D}}\cap\partial\Omega^{\mathrm{E}}$. On
$\partial\Omega^{\mathrm{E}}$ we define a part of the boundary related to the
exposed surface, denoted as $\Gamma$ (see Figure 2.1).

$\partial\Omega^{\mathrm{D}}$$\partial\Omega^{\mathrm{E}}$ Figure 2.1:
Schematic representation of dermis and epidermis layers separated by the
interface $\Sigma$ and an exposed epidermic surface is denoted by $\Gamma$.

This heterogeneity of the skin (corresponding to the presence of an
homogeneous _dermis_ and an homogeneous _epidermis_ , occupied by
$\Omega^{\mathrm{D}}$ and $\Omega^{\mathrm{E}}$, respectively) is accounted
for by the spatially-dependent (but locally constant) Lamé coefficients

$\lambda=\begin{cases}\lambda^{\mathrm{D}}&\text{in $\Omega^{\mathrm{D}}$},\\\
\lambda^{\mathrm{E}}&\text{in
$\Omega^{\mathrm{E}}$},\end{cases}\quad\mu=\begin{cases}\mu^{\mathrm{D}}&\text{in
$\Omega^{\mathrm{D}}$},\\\ \mu^{\mathrm{E}}&\text{in
$\Omega^{\mathrm{E}}$}.\end{cases}$

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 $\boldsymbol{u}^{*}$ and the pressure
$p^{*}=-\lambda^{*}\operatorname*{div}\boldsymbol{u}^{*}$, with
$*\in\\{$D,E$\\}$ (although the interpretation of $p$ as pressure holds
exclusively in the incompressibility limit where $\lambda^{*}/\mu^{*}$ takes
on quite large values). For all $0<t<T$, and given a body force
$\boldsymbol{f}$ (whose specific form depends on the vector $\boldsymbol{w}$;
see equation (2.5) below) applied on each domain, we resolve the displacement
/ pressure pair
$(\boldsymbol{u}^{*}(t),p^{*}(t)):\Omega^{*}\to[\mathbb{R}^{d}\times\mathbb{R}]$
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 $\Sigma$). Then the
global linear Navier-Lamé problem adopts the form

$\displaystyle-\operatorname*{\mathbf{div}}\bigl{(}2\mu^{*}\boldsymbol{\varepsilon}(\boldsymbol{u}^{*})-p^{*}\mathbf{I}\bigr{)}$
$\displaystyle=\boldsymbol{f}(\boldsymbol{w}^{*})$ $\displaystyle\text{in
$\Omega^{*}\times(0,T)$},$ $\displaystyle
p^{*}+\lambda^{*}\operatorname*{div}\boldsymbol{u}^{*}$ $\displaystyle=0$
$\displaystyle\text{in $\Omega^{*}\times(0,T)$},$
$\displaystyle\boldsymbol{\sigma}^{\mathrm{E}}\boldsymbol{\nu}+\alpha^{\mathrm{E}}\boldsymbol{u}^{\mathrm{E}}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\text{on $\Gamma\times(0,T)$},$
(2.1) $\displaystyle\boldsymbol{u}^{*}$ $\displaystyle=\boldsymbol{0}$
$\displaystyle\text{on $\partial\Omega\setminus\Gamma\times(0,T)$},$
$\displaystyle\boldsymbol{\sigma}^{\mathrm{D}}\boldsymbol{\nu}=\boldsymbol{\sigma}^{\mathrm{E}}\boldsymbol{\nu},$
$\displaystyle\quad\boldsymbol{u}^{\mathrm{D}}=\boldsymbol{u}^{\mathrm{E}}$
$\displaystyle\text{on $\Sigma\times(0,T)$},$

where
$\boldsymbol{\varepsilon}(\boldsymbol{u}^{*})=\frac{1}{2}(\nabla\boldsymbol{u}^{*}+\nabla\boldsymbol{u}^{*T})$
is the tensor of infinitesimal strains (symmetric gradient of displacement),
$\boldsymbol{\sigma}^{*}=2\mu^{*}\boldsymbol{\varepsilon}(\boldsymbol{u}^{*})-p^{*}\mathbf{I}$
is the Cauchy stress tensor, and $\alpha^{*}>0$ is the Robin coefficient. The
extension to non-homogeneous displacement conditions
($\boldsymbol{u}|_{\partial\Omega}\neq 0$) can be handled in a standard way by
adding suitable liftings on the data. Below, we use the asterisk superscript
$*\in\\{$D,E$\\}$ to denote quantities related to each subdomain, and we adopt
the convention that the unit normal vector $\boldsymbol{\nu}$ on the interface
is pointing from $\Omega^{\mathrm{D}}$ to $\Omega^{\mathrm{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 $\boldsymbol{w}=(w_{1},\ldots,w_{m})$. 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 $\Gamma$
and continuity of fluxes and densities at the interface) leads to the system

$\displaystyle\partial_{t}\boldsymbol{w}^{*}+(\partial_{t}\boldsymbol{u}^{*}\cdot\nabla)\boldsymbol{w}^{*}-\operatorname*{\mathbf{div}}(\mathbf{M}^{*}\nabla\boldsymbol{w}^{*})$
$\displaystyle=\mathbf{G}(\boldsymbol{w}^{*})+\boldsymbol{g}(\boldsymbol{u}^{*})$
$\displaystyle\text{in $\Omega^{*}\times(0,T)$},$
$\displaystyle\mathbf{M}^{*}\nabla\boldsymbol{w}^{*}\boldsymbol{\nu}$
$\displaystyle=\boldsymbol{0}$ $\displaystyle\text{on
$\partial\Omega\times(0,T)$},$ (2.2)
$\displaystyle\mathbf{M}^{\mathrm{D}}\nabla\boldsymbol{w}^{\mathrm{D}}\boldsymbol{\nu}=\mathbf{M}^{\mathrm{E}}\nabla\boldsymbol{w}^{\mathrm{E}}\boldsymbol{\nu},$
$\displaystyle\quad\boldsymbol{w}^{\mathrm{D}}=\boldsymbol{w}^{\mathrm{E}}$
$\displaystyle\text{on $\Sigma\times(0,T)$}.$

Here $\mathbf{M}^{*}\in\mathbb{R}^{m\times m}$ is a (not necessarily diagonal,
but positive definite) tensor of self- and cross-diffusion rates,
$\mathbf{G}:\Omega\to\mathbb{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
$\Omega^{\mathrm{D}},\Omega^{\mathrm{E}}$; and the term $\boldsymbol{g}$
encodes an additional coupling contribution (detailed in equation (2.5)
below). The system is endowed with adequate initial conditions for
$\boldsymbol{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 activator-inhibitor model with quadratic reaction terms due to
Gierer-Meinhardt (cf. [24]):

$\mathbf{G}(\boldsymbol{w}^{*})=\begin{pmatrix}\rho_{2}(\rho_{0}+\rho_{1}\frac{(w_{1}^{*})^{2}}{w^{*}_{2}})-\rho_{3}w^{*}_{1}\\\
\rho_{4}(w^{*}_{1})^{2}-\rho_{5}w^{*}_{2}\end{pmatrix},$ (2.3)

where $w^{*}_{1}$ is the short-range autocatalytic substance and $w^{*}_{2}$
is the long-range antagonist; and the $\rho_{i}$’s are positive model
parameters.

Second, we introduce a more general model of skin appendages development were
$w^{*}_{1}$ is the cell density, $w^{*}_{2}$ is the density of the
extracellular matrix, and $w^{*}_{3},w_{4}^{*}$ are the concentrations of two
morphogen signalling chemical compounds. All quantities are assumed to
diffuse, so none of the diagonal terms in $\mathbf{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 $\mathbf{M}^{\mathrm{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:

$\mathbf{G}(\boldsymbol{w}^{\mathrm{D}})=\begin{pmatrix}r_{1}w_{1}^{\mathrm{D}}(r_{0}-w_{1}^{\mathrm{D}})\\\
0\\\ r_{2}w_{1}^{\mathrm{D}}-r_{3}w_{3}^{\mathrm{D}}\\\
-r_{4}w_{1}^{\mathrm{D}}w_{4}^{\mathrm{D}}\end{pmatrix},\quad\mathbf{G}(\boldsymbol{w}^{\mathrm{E}})=\begin{pmatrix}0\\\
0\\\ -r_{5}w_{1}^{\mathrm{E}}w_{3}^{\mathrm{E}}\\\
r_{6}w_{1}^{\mathrm{E}}-r_{7}w_{4}^{\mathrm{E}}\end{pmatrix},$ (2.4)

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]).

#### Coupling mechanisms.

Extensive experimental evidence indicates that both biochemical and mechanical
processes work in concert during pattern formation (see for instance [1, 2,
25, 30, 35, 53]). 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 $\boldsymbol{g}(\boldsymbol{u}^{*})$ on the right-hand
side in (2.2), carrying local information about the medium dilation.
Therefore, we set

$\boldsymbol{f}(\boldsymbol{w}^{*})=c_{f}^{*}\sum_{i=1}^{m}\nabla
w^{*}_{i},\qquad\boldsymbol{g}(\boldsymbol{u}^{*})=c_{g}^{*}\operatorname*{div}\boldsymbol{u}^{*}\begin{pmatrix}1\\\
\vdots\\\ 1\end{pmatrix},$ (2.5)

where $c_{f}^{*},c_{g}^{*}$ are model parameters. Note that the form of
$\boldsymbol{f}$ is equivalent to considering a stress component depending on
the concentrations of some $w_{i}^{*}$, that is
$\boldsymbol{\sigma}^{*}=\boldsymbol{\sigma}^{*}_{\mathrm{eff}}+w^{*}_{i}\mathbf{I}$,
where $\boldsymbol{\sigma}^{*}_{\mathrm{eff}}$ is an effective elastic stress.
Similar coupling conditions can be found in [22, 31, 49].

###  2.1 Weak formulation for the multidomain problem

Multiplying (2.1)-(2.2) by suitable test functions, integrating by parts over
each region $\Omega^{*}$, and exploiting the boundary and interface
conditions, we arrive at the following primal-mixed variational problem
associated to the elasticity-ADR problem: For $t>0$, find
$(\boldsymbol{u}^{*}(t),p^{*}(t),\boldsymbol{w}^{*}(t))\in\boldsymbol{\mathcal{V}}^{*}\times\mathcal{Q}^{*}\times\boldsymbol{\mathcal{Z}}^{*}$
such that

$\displaystyle
a(\boldsymbol{u}^{*}(t),\boldsymbol{v})+b(p^{*}(t),\boldsymbol{v})+\alpha^{*}\langle\boldsymbol{u}^{*}(t),\boldsymbol{v}\rangle_{\Gamma}$
$\displaystyle=F(\boldsymbol{w}^{*};\boldsymbol{v})$
$\displaystyle\forall\boldsymbol{v}\in\boldsymbol{\mathcal{V}}^{*},$
$\displaystyle b\bigl{(}q,\boldsymbol{u}^{*}(t)\bigl{)}-c(p^{*}(t),q)$
$\displaystyle=0$ $\displaystyle\forall q\in\mathcal{Q}^{*},$ (2.6)
$\displaystyle\bigl{(}\partial_{t}\boldsymbol{w}^{*}(t),\boldsymbol{z}\bigr{)}_{\Omega^{*}}+C(\boldsymbol{u}^{*}(t),\boldsymbol{w}^{*}(t);\boldsymbol{z})+d(\boldsymbol{w}^{*}(t),\boldsymbol{z})$
$\displaystyle=G(\boldsymbol{w}^{*};\boldsymbol{z})$
$\displaystyle\forall\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}^{*},$

where $\alpha^{\mathrm{E}}$ is the coefficient associated to the Robin
boundary condition at the surface, and the specific boundary treatment implies
that adequate functional spaces are

$\boldsymbol{\mathcal{V}}^{*}:=H^{1}_{0}(\Omega^{*})^{d}=\\{\boldsymbol{v}\in
H^{1}(\Omega^{*})^{d}:\boldsymbol{v}|_{\partial\Omega\setminus\Gamma}=\boldsymbol{0}\\},\qquad\mathcal{Q}^{*}=L^{2}(\Omega^{*}),\qquad\boldsymbol{\mathcal{Z}}^{*}=H^{1}(\Omega^{*})^{m}.$

Linear and nonlinear forms (lowercase and uppercase, respectively) involved in
(2.6) are defined as

$\displaystyle a(\boldsymbol{u},\boldsymbol{v})$
$\displaystyle:=2\mu^{*}(\boldsymbol{\varepsilon}(\boldsymbol{u}),\boldsymbol{\varepsilon}(\boldsymbol{v}))_{\Omega^{*}},\quad
b(q,\boldsymbol{v}):=-(q,\operatorname*{div}\boldsymbol{v})_{\Omega^{*}},\quad
c(p,q):=\frac{1}{\lambda^{*}}(p,q)_{\Omega^{*}},$ $\displaystyle
C(\boldsymbol{u}(t),\boldsymbol{w};\boldsymbol{z})$
$\displaystyle:=\bigl{(}[\partial_{t}\boldsymbol{u}\cdot\nabla]\boldsymbol{w},\boldsymbol{z}\bigr{)}_{\Omega^{*}},\quad
d(\boldsymbol{w},\boldsymbol{z}):=(\mathbf{M}^{*}\nabla\boldsymbol{w},\nabla\boldsymbol{z})_{\Omega^{*}},$
$\displaystyle F(\boldsymbol{w};\boldsymbol{v})$
$\displaystyle:=\bigl{(}\boldsymbol{f}(\boldsymbol{w}),\boldsymbol{v}\bigr{)}_{\Omega^{*}},\quad
G(\boldsymbol{u},\boldsymbol{w};\boldsymbol{z}):=\bigl{(}\mathbf{G}(\boldsymbol{w})+g(\boldsymbol{u}),\boldsymbol{z}\bigr{)}_{\Omega^{*}}.$

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 $\Omega$. 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

$\mathcal{R}(\boldsymbol{w}^{\mathrm{D}}):=\mathbf{M}^{\mathrm{D}}\nabla\boldsymbol{w}^{\mathrm{D}}\boldsymbol{\nu}+K^{\mathrm{D}}\boldsymbol{w}^{\mathrm{D}}=\mathbf{M}^{\mathrm{E}}\nabla\boldsymbol{w}^{\mathrm{E}}\boldsymbol{\nu}+K^{\mathrm{E}}\boldsymbol{w}^{\mathrm{E}}=:\mathcal{R}(\boldsymbol{w}^{\mathrm{E}})\qquad\text{on
$\Sigma\times(0,T)$},$ (2.7)

where $\mathcal{R}(\cdot)$ denotes the Robin transmission operator, and
$K^{\mathrm{D}},K^{\mathrm{E}}$ are nonnegative acceleration constants
satisfying $K^{\mathrm{D}}+K^{\mathrm{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

$\mathcal{S}(\boldsymbol{u}^{\mathrm{D}},p^{\mathrm{D}}):=\boldsymbol{\sigma}^{\mathrm{D}}\boldsymbol{\nu}+J^{\mathrm{D}}\boldsymbol{u}^{\mathrm{D}}=\boldsymbol{\sigma}^{\mathrm{E}}\boldsymbol{\nu}+J^{\mathrm{E}}\boldsymbol{u}^{\mathrm{E}}=:\mathcal{S}(\boldsymbol{u}^{\mathrm{E}},p^{\mathrm{E}})\qquad\text{on
$\Sigma\times(0,T)$},$ (2.8)

where $\mathcal{S}(\cdot,\cdot)$ 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 $\ell=0,\ldots$, find
$(\boldsymbol{u}^{*,\ell+1},p^{*,\ell+1},\boldsymbol{w}^{*,\ell+1})\in\boldsymbol{\mathcal{V}}^{*}\times\mathcal{Q}^{*}\times\boldsymbol{\mathcal{Z}}^{*}$
such that

$\displaystyle
a(\boldsymbol{u}^{*,\ell+1},\boldsymbol{v})+b(p^{*,\ell+1},\boldsymbol{v})+L^{*}\langle\boldsymbol{u}^{*,\ell+1},\boldsymbol{v}\rangle_{\Sigma\cup\Gamma}$
$\displaystyle=F(\boldsymbol{w}^{*,\ell+1};\boldsymbol{v})+\langle\mathcal{S}^{{\dagger},\ell},\boldsymbol{v}\rangle_{\Sigma}$
$\displaystyle\forall\boldsymbol{v}\in\boldsymbol{\mathcal{V}}^{*},$
$\displaystyle b\bigl{(}q,\boldsymbol{u}^{*,\ell+1}\bigl{)}-c(p^{*,\ell+1},q)$
$\displaystyle=0$ $\displaystyle\forall q\in\mathcal{Q}^{*},$ (2.9)
$\displaystyle\bigl{(}\partial_{t}\boldsymbol{w}^{*,\ell+1},\boldsymbol{z}\bigr{)}_{\Omega^{*}}+C(\boldsymbol{u}^{*,\ell+1},\boldsymbol{w}^{*,\ell+1};\boldsymbol{z})\qquad$
$\displaystyle+d(\boldsymbol{w}^{*,\ell+1},\boldsymbol{z})+K^{*}\langle\boldsymbol{w}^{*,\ell+1},\boldsymbol{z}\rangle_{\Sigma}$
$\displaystyle=G(\boldsymbol{w}^{*,\ell+1};\boldsymbol{z})+\langle\mathcal{R}(\boldsymbol{w}^{{\dagger},\ell}),\boldsymbol{z}\rangle_{\Sigma}$
$\displaystyle\forall\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}^{*},$

where

$L^{*}=\begin{cases}\alpha^{\mathrm{E}}&\text{if $*=$\,E and on $\Gamma$},\\\
J^{*}&\text{ on $\Sigma$},\\\ 0&\text{otherwise},\end{cases}$

and the superscript ${\dagger}\in\\{$D,E$\\}$ denotes quantities in the
subdomain opposite to $*$. Here we have dropped the explicit time-dependence
whenever clear from the context.

##  3 A finite element method and its splitting scheme

In this section we detail a Rothe-type strategy employed to numerically solve
the set of governing equations (2.1)-(2.2). We first apply a discretisation in
space written as a primal-mixed Galerkin formulation, which we later combine
with fixed-point schemes and a Schwarz iterative algorithm decoupling the
problem into four main blocks corresponding to each side of the interface, and
to the elasticity and ADR equations. Then we discuss the time discretisation
based on an adaptive scheme and operator splitting.

###  3.1 Nonlinear, partitioned Galerkin formulation

Let us denote by $\mathcal{T}_{h}^{*}$ a shape-regular simplicial partition of
$\Omega^{*}$ by triangles (for $d=2$, or tetrahedra for $d=3$) $K$ of diameter
$h_{K}$ and define the global meshsize as
$h:=\max\\{h_{K}:K\in\mathcal{T}^{*}_{h}\\}$. For clarity of the presentation,
we assume that the two meshes (associated with each subdomain) match at the
interface. We will seek the approximation
$(\boldsymbol{u}_{h}^{*},p_{h}^{*},\boldsymbol{w}_{h}^{*})$ of
$(\boldsymbol{u}^{*},p^{*},\boldsymbol{w}^{*})$, in the finite dimensional
product space

$[\boldsymbol{\mathcal{V}}^{*}_{h}\times\mathcal{Q}^{*}_{h}\times\boldsymbol{\mathcal{Z}}^{*}_{h}]\subset[\boldsymbol{\mathcal{V}}^{*}\times\mathcal{Q}^{*}\times\boldsymbol{\mathcal{Z}}^{*}],$
(3.1)

to be made precise below. Then, the semidiscrete Galerkin scheme associated to
the Schwarz alternating algorithm (2.9) and a fixed-point decoupling of
elasticity and ADR subproblems can be recast in the following partitioned
form, iterated until convergence: For a given $t>0$, and given
$(\boldsymbol{u}_{h}^{*,\ell},p_{h}^{*,\ell},\boldsymbol{w}_{h}^{*,\ell})\in\boldsymbol{\mathcal{V}}_{h}^{*}\times\mathcal{Q}_{h}^{*}\times\boldsymbol{\mathcal{Z}}_{h}^{*}$,
find the next iteration
$(\boldsymbol{u}_{h}^{*,\ell+1},p_{h}^{*,\ell+1},\boldsymbol{w}_{h}^{*,\ell+1})$
from the sequential system

$\displaystyle\begin{cases}&\bigl{(}\partial_{t}\boldsymbol{w}_{h}^{\mathrm{D},\ell+1},\boldsymbol{z}_{h}^{\mathrm{D}}\bigr{)}_{\Omega^{\mathrm{D}}}+C(\boldsymbol{u}_{h}^{\mathrm{D},\ell},\boldsymbol{w}_{h}^{\mathrm{D},\ell+1};\boldsymbol{z}_{h}^{\mathrm{D}})+d(\boldsymbol{w}_{h}^{\mathrm{D},\ell+1},\boldsymbol{z}_{h})\\\
&\qquad+K^{\mathrm{D}}\langle\boldsymbol{w}_{h}^{\mathrm{D},\ell+1},\boldsymbol{z}_{h}^{\mathrm{D}}\rangle_{\Sigma}-G(\boldsymbol{w}_{h}^{\mathrm{D},\ell+1};\boldsymbol{z}_{h}^{\mathrm{D}})=\langle\mathcal{R}_{h}^{\mathrm{E},\ell},\boldsymbol{z}_{h}^{\mathrm{D}}\rangle_{\Sigma}\quad\forall\boldsymbol{z}_{h}^{\mathrm{D}}\in\boldsymbol{\mathcal{Z}}_{h}^{\mathrm{D}},\end{cases}$
(3.2)
$\displaystyle\begin{cases}a(\boldsymbol{u}_{h}^{\mathrm{D},\ell+1},\boldsymbol{v}_{h}^{\mathrm{D}})+b(p_{h}^{\mathrm{D},\ell+1},\boldsymbol{v}_{h}^{\mathrm{D}})+L^{\mathrm{D}}\langle\boldsymbol{u}_{h}^{\mathrm{D},\ell+1},\boldsymbol{v}_{h}^{\mathrm{D}}\rangle_{\Sigma\cup\Gamma}\\!\\!\\!&=F(\boldsymbol{w}_{h}^{\mathrm{D},\ell};\boldsymbol{v}_{h}^{\mathrm{D}})+\langle\mathcal{S}_{h}^{\mathrm{E},\ell},\boldsymbol{v}_{h}^{\mathrm{D}}\rangle_{\Sigma}\
\forall\boldsymbol{v}_{h}^{\mathrm{D}}\in\boldsymbol{\mathcal{V}}_{h}^{\mathrm{D}},\\\
b\bigl{(}q_{h}^{\mathrm{D}},\boldsymbol{u}_{h}^{\mathrm{D},\ell+1}\bigl{)}-\,c(p_{h}^{\mathrm{D},\ell+1},q_{h}^{\mathrm{D}})&=0\quad\forall
q_{h}^{\mathrm{D}}\in\mathcal{Q}_{h}^{\mathrm{D}},\end{cases}$ (3.3)
$\displaystyle\begin{cases}&\bigl{(}\partial_{t}\boldsymbol{w}_{h}^{\mathrm{E},\ell+1},\boldsymbol{z}_{h}\mathrm{E}\bigr{)}_{\Omega^{\mathrm{E}}}+C(\boldsymbol{u}_{h}^{\mathrm{E},\ell},\boldsymbol{w}_{h}^{\mathrm{E},\ell+1};\boldsymbol{z}_{h}^{\mathrm{E}})+d(\boldsymbol{w}_{h}^{\mathrm{E},\ell+1},\boldsymbol{z}_{h}^{\mathrm{E}})\\\
&\qquad+K^{\mathrm{E}}\langle\boldsymbol{w}_{h}^{\mathrm{E},\ell+1},\boldsymbol{z}_{h}^{\mathrm{E}}\rangle_{\Sigma}-G(\boldsymbol{w}_{h}^{\mathrm{E},\ell+1};\boldsymbol{z}_{h}^{\mathrm{E}})=\langle\mathcal{R}_{h}^{\mathrm{D},\ell+1},\boldsymbol{z}_{h}^{\mathrm{E}}\rangle_{\Sigma}\quad\forall\boldsymbol{z}_{h}^{\mathrm{E}}\in\boldsymbol{\mathcal{Z}}_{h}^{\mathrm{E}},\end{cases}$
(3.4)
$\displaystyle\begin{cases}a(\boldsymbol{u}_{h}^{\mathrm{E},\ell+1},\boldsymbol{v}_{h}^{\mathrm{E}})+b(p_{h}^{\mathrm{E},\ell+1},\boldsymbol{v}_{h}^{\mathrm{E}})+L^{\mathrm{E}}\langle\boldsymbol{u}_{h}^{\mathrm{E},\ell+1},\boldsymbol{v}_{h}^{\mathrm{E}}\rangle_{\Sigma\cup\Gamma}\\!\\!\\!\\!&=F(\boldsymbol{w}_{h}^{\mathrm{E},\ell+1};\boldsymbol{v}_{h}^{\mathrm{E}})+\langle\mathcal{S}_{h}^{\mathrm{D},\ell+1},\boldsymbol{v}_{h}^{\mathrm{E}}\rangle_{\Sigma}\
\forall\boldsymbol{v}_{h}^{\mathrm{E}}\in\boldsymbol{\mathcal{V}}_{h}^{\mathrm{E}},\\\
b\bigl{(}q_{h}^{\mathrm{E}},\boldsymbol{u}_{h}^{\mathrm{E},\ell+1}\bigl{)}-\,c(p_{h}^{\mathrm{E},\ell+1},q_{h}^{\mathrm{E}})&=0\quad\forall
q_{h}^{\mathrm{E}}\in\mathcal{Q}_{h}^{\mathrm{E}}.\end{cases}$ (3.5)

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

$\boldsymbol{\mathcal{V}}^{*}_{h}:=\\{\boldsymbol{v}\in\boldsymbol{\mathcal{V}}^{*}:\boldsymbol{v}|_{K}\in[\mathbb{P}_{1}(K)+\mathbb{B}(K)]^{d},\
K\in\mathcal{T}_{h}^{*}\\},\quad\mathcal{Q}_{h}^{*}:=\\{q\in\mathcal{Q}^{*}:q|_{K}\in\mathbb{P}_{1}(K),\
K\in\mathcal{T}_{h}^{*}\\},$

where $\mathbb{B}(K)$ 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
$\boldsymbol{\mathcal{Z}}_{h}^{*}:=\\{\boldsymbol{z}\in\boldsymbol{\mathcal{Z}}^{*}:\boldsymbol{z}|_{K}\in\mathbb{P}_{1}(K)^{m},\
K\in\mathcal{T}_{h}^{*}\\}.$

The discretisation (3.2)-(3.5) leads to the following system of equations

$\displaystyle\mathbfsf{M}^{*}\dot{\mathbf{w}}^{*,l+1}+\mathbfsf{C}_{u}^{*}\mathbf{w}^{*,l+1}+\mathbfsf{D}\mathbf{w}^{*,l+1}-\mathbf{G}^{*}(\mathbf{w}^{*,l+1})$
$\displaystyle=\boldsymbol{0},$ (3.6)
$\displaystyle\begin{cases}\mathbfsf{A}^{*}\mathbf{u}^{*,l+1}+(\mathbfsf{B}^{*})^{T}\mathbf{p}^{*,l+1}-\mathbf{F}^{*}(\mathbf{w}^{*,\ddagger})&=\boldsymbol{0},\\\
\mathbfsf{B}^{*}\mathbf{u}^{*,l+1}-\mathbfsf{C}\mathbf{p}^{*,l+1}&=\boldsymbol{0},\end{cases}$
(3.7)

where the superscript $\ddagger$ 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 $\mathbfsf{C}_{u}^{*}$
emphasises that the matrix is directly dependent on displacements.

###  3.2 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
$\frac{d}{dt}v_{h}=L_{h}(v_{h},\varphi_{h}(t))$, we recall that the main step
in the classical RK23 scheme reads

$\displaystyle v_{h}^{(i)}$
$\displaystyle=v_{h}^{n}+\sum_{j=1}^{s}a_{ij}\Delta
t^{n}L_{h}(v_{h}^{(j)},\varphi_{h}(t^{n}+c_{j}\Delta t^{n}))\quad i=1,2,3,$
$\displaystyle v_{h}^{n+1}$ $\displaystyle=v_{h}^{n}+\sum_{j=1}^{s}b_{j}\Delta
t^{n}L_{h}(v_{h}^{(j)},\varphi_{h}(t^{n}+c_{j}\Delta t^{n})),$
$\displaystyle\widehat{v}_{h}^{n+1}$
$\displaystyle=\widehat{v}_{h}^{n}+\sum_{j=1}^{s}\widehat{b}_{j}\Delta
t^{n}L_{h}(v_{h}^{(j)},\varphi_{h}(t^{n}+c_{j}\Delta t^{n})),$

with $s=3$ stages, and where the two approximations $v_{h},\widehat{v}_{h}$
are of order $\mathcal{O}(p)$ and $\mathcal{O}(\widehat{p})$, respectively. As
in classical embedded RK methods, here the coefficients $\\{a_{ij},c_{j}\\}$
coincide in both approximations; however the coefficients $b_{j}$ differ in
order to attain the expected convergence orders (in general one chooses
$\widehat{p}=p\pm 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.

0 | 0 | 0 | 0
---|---|---|---
$2\gamma$ | $\gamma$ | $\gamma$ | 0
1 | $1-b_{2}-\gamma$ | $b_{2}$ | $\gamma$
| $1-b_{2}-\gamma$ | $b_{2}$ |  $\gamma$
| $1-\widehat{b}_{2}-\widehat{b}_{3}$ | $\widehat{b}_{2}$ |  $\widehat{b}_{3}$
Table 3.1: Butcher’s tableau for the TR-BDF2 method.

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}\neq 0$) [27]. The scheme is
characterised by the following Butcher’s tableau reported in Table 3.1, where
$b_{2}=\frac{1-2\gamma}{4\gamma}$,
$\widehat{b}_{2}=\frac{1}{12\gamma(1-2\gamma)}$,
$\widehat{b}_{3}=\frac{2-6\gamma}{6(1-2\gamma)}$ and
$\gamma=\frac{2-\sqrt{2}}{2}$.

As established in [29], the A-stability, the stiffly accuracy ($a_{sj}=b_{j}$,
for $j=1,\ldots,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.

Based on Table 3.1, and considering
$L_{h}=\mathbfsf{C}_{u}^{*}\mathbf{w}^{*}+\mathbfsf{D}\mathbf{w}^{*}-\mathbf{G}^{*}(\mathbf{w}^{*})$
(where we have removed the Schwarz index whenever clear from the context), we
can recast the intermediate stages as

$\displaystyle\mathbf{w}^{*,(1)}$ $\displaystyle=\mathbf{w}^{*,n},$
$\displaystyle\mathbfsf{M}^{*}\mathbf{w}^{*,(2)}$ $\displaystyle=\Delta
t^{n}\gamma\bigl{[}-\mathbfsf{C}_{u}^{*}\mathbf{w}^{*,(2)}-\mathbfsf{D}\mathbf{w}^{*,(2)}+\mathbf{G}^{*}(\mathbf{w}^{*,(2)})\bigr{]}$
$\displaystyle\qquad+\Delta
t^{n}\gamma\bigl{[}-\mathbfsf{C}_{u}^{*}\mathbf{w}^{*,(1)}-\mathbfsf{D}\mathbf{w}^{*,(1)}+\mathbf{G}^{*}(\mathbf{w}^{*,(1)})\bigr{]}+\mathbfsf{M}^{*}\mathbf{w}^{*,n},$
$\displaystyle\mathbfsf{M}^{*}\mathbf{w}^{*,(3)}$ $\displaystyle=\Delta
t^{n}\gamma\bigl{[}-\mathbfsf{C}_{u}^{*}\mathbf{w}^{*,(3)}-\mathbfsf{D}\mathbf{w}^{*,(3)}+\mathbf{G}^{*}(\mathbf{w}^{*,(3)})\bigr{]}$
$\displaystyle\qquad+\Delta
t^{n}b_{2}\bigl{[}-\mathbfsf{C}_{u}^{*}\mathbf{w}^{*,(2)}-\mathbfsf{D}\mathbf{w}^{*,(2)}+\mathbf{G}^{*}(\mathbf{w}^{*,(2)})\bigr{]}$
$\displaystyle\qquad+\Delta
t^{n}(1-b_{2}-\gamma)\bigl{[}-\mathbfsf{C}_{u}^{*}\mathbf{w}^{*,(1)}-\mathbfsf{D}\mathbf{w}^{*,(1)}+\mathbf{G}^{*}(\mathbf{w}^{*,(1)})\bigr{]}+\mathbfsf{M}^{*}\mathbf{w}^{*,n}.$

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_{z,h}^{(i)}:=v_{h}^{(i)}-v_{h}^{n}$ and rewrite the
problem accordingly. Doing so results in decreasing round-off errors and
serves in reducing the cost to compute $\widehat{v}_{h}^{n+1}$. We follow this
approach and modify the superscripts $(1)\rightarrow n$, $(2)\rightarrow
n+\gamma$ and $(3)\rightarrow new$, to obtain the following full
discretisation of the ODE system:

$\displaystyle\mathbf{w}^{*,n+\gamma}$
$\displaystyle=\gamma\mathbf{w}_{z}^{*,n+\gamma}+\gamma\mathbf{w}_{z}^{*,n}+\mathbf{w}^{*,n},$
$\displaystyle\mathbf{w}^{*,new}$
$\displaystyle=\gamma\mathbf{w}_{z}^{*,new}+b_{2}\mathbf{w}_{z}^{*,n+\gamma}+(1-b_{2}-\gamma)\mathbf{w}_{z}^{*,n}+\mathbf{w}^{*,n},$

with

$\displaystyle\mathbf{w}_{z}^{*,n+\gamma}$
$\displaystyle=(\mathbfsf{M}^{*})^{-1}\Delta
t^{n}\left[-\mathbfsf{C}_{u}^{*}\mathbf{w}^{*,n+\gamma}-\mathbfsf{D}\mathbf{w}^{*,n+\gamma}+\mathbf{G}^{*}(\mathbf{w}^{*,n+\gamma})\right],$
$\displaystyle\mathbf{w}_{z}^{*,new}$
$\displaystyle=(\mathbfsf{M}^{*})^{-1}\Delta
t^{n}\left[-\mathbfsf{C}_{u}^{*}\mathbf{w}^{*,new}-\mathbfsf{D}\mathbf{w}^{*,new}+\mathbf{G}^{*}(\mathbf{w}^{*,new})\right].$

This nonlinear system is solved by means of Newton steps as follows

$\displaystyle\mathbf{w}_{z,k+1}^{*}$
$\displaystyle=\mathbf{w}_{z,k}^{*}+\delta_{k+1}\mathbf{w}_{z}^{*},$
$\displaystyle\biggl{(}\frac{1}{\Delta
t^{n}}\mathbfsf{M}^{*}-\gamma\mathbfsf{J}^{*}\biggr{)}^{-1}\delta_{k+1}\mathbf{w}_{z}^{*}$
$\displaystyle=\left[-\mathbfsf{C}_{u}^{*}\mathbf{w}_{k}^{*}-\mathbfsf{D}\mathbf{w}_{k}^{*}+\mathbf{G}^{*}(\mathbf{w}_{k}^{*})\right]-\frac{1}{\Delta
t^{n}}\mathbfsf{M}\mathbf{w}_{z,k}^{*},$

where

$\mathbfsf{J}^{*}=-\mathbfsf{D}^{*}+\frac{\partial\mathbf{G}^{*}(\mathbf{w}_{k}^{*})}{\partial\mathbf{w}^{*}}-\mathbfsf{C}_{u}^{*}=-\mathbfsf{D}^{*}+\mathbfsf{J}_{\mathbf{G}^{*}}-\mathbfsf{C}_{u}^{*}.$

We finally point out that an additional reduction of computational cost occurs
due to the fact that the local error
$\mathrm{err}=\widehat{v}_{h}^{n+1}-v_{h}^{n+1}$ can be accessed by the simple
relation

$\mathrm{err}=(\widehat{b}_{1}-b_{1})\mathbf{w}_{z}^{n}+(\widehat{b}_{2}-b_{2})\mathbf{w}_{z}^{n+\gamma}+(\widehat{b}_{3}-b_{3})\mathbf{w}_{z}^{new},$

avoiding the explicit storage and computation of $\widehat{v}_{h}^{n+1}$.

###  3.3 Algorithmic details

Overall solution algorithm. The global numerical scheme to solve the
mechanochemical model is summarised in Alg. 1. After assigning suitable
initial conditions $\mathbf{u}^{*,0}=\boldsymbol{0}$,
$\mathbf{p}^{*,0}=\boldsymbol{0}$, $\mathbf{w}^{*,0}=\mathbf{w}_{0}^{*}$,
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 $\mathbf{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.

Algorithm 1 Overall solution algorithm for the mechanochemical problem

1:Initialise the solution vectors $\mathbf{u}^{*,0}$, $\mathbf{p}^{*,0}$,
$\mathbf{w}^{*,0}$ and define
$\mathbf{w}\triangleq\left[\mathbf{w}^{D},\mathbf{w}^{E}\right]^{T}$,
$\mathbf{w}_{z}\triangleq\left[\mathbf{w}_{z}^{D},\mathbf{w}_{z}^{E}\right]^{T}$

2:while $t\leq T_{f}$ do

3: $\Delta t^{\min}=\epsilon|t|$, $\Delta t^{\mathrm{old}}=\Delta t$, $\Delta
t=\min(\Delta t^{\max},\max(\Delta t^{\min},\Delta t^{\mathrm{old}}))$

4: if $\Delta t\neq\Delta t^{\mathrm{old}}$ then

5: $\mathbf{w}_{z}^{n}=\frac{\Delta t}{\Delta
t^{\mathrm{old}}}\mathbf{w}_{z}^{n}$

6: needNewMJ $\leftarrow$ true

7: end if

8: Compute
$t_{\mathrm{new}},\mathbf{w}^{\mathrm{new}},\mathbf{w}_{z}^{\mathrm{new}}$
using TR-BDF2($\mathbf{w}^{n}$,$\mathbf{w}_{z}^{n}$,$\Delta
t$,$\boldsymbol{\beta}_{\Delta t}$) $\triangleright$ See Algorithm 2

9: for $*\in\\{D,E\\}$ do procedure ElasticitySolver associated to (3.3) and
(3.5)

10: $\mathbf{F}^{*,n+1}=\mathbf{F}(\mathbf{w}^{*,new})$,
$\mathbf{v}^{*,n+1}=-\frac{\mathbf{u}^{*,n}}{\Delta t}$

11: $\left[\begin{array}[]{c}\mathbf{u}^{*,n+1}\\\
\mathbf{p}^{*,n+1}\end{array}\right]=\left[\begin{array}[]{cc}\mathbfsf{A}^{*}&(\mathbfsf{B}^{*})^{T}\\\
\mathbfsf{B}^{*}&-\mathbfsf{C}^{*}\\\
\end{array}\right]^{-1}\left[\begin{array}[]{c}\mathbf{F}^{*,n+1}\\\
\boldsymbol{0}\end{array}\right]$

12: $\mathbf{v}^{*,n+1}=\mathbf{v}^{*,n+1}+\frac{\mathbf{u}^{*,n+1}}{\Delta
t}$

13: update Robin transmission conditions $\mathcal{S}_{h}^{*}$ from (2.8)

14: end for

15: set $n=n+1$, $t_{n+1}=t_{\mathrm{new}}$,
$\mathbf{w}^{n+1}=\mathbf{w}^{\mathrm{new}}$,
$\mathbf{w}_{z}^{n+1}=\mathbf{w}_{z}^{\mathrm{new}}$, and
$|\mathbf{w}^{n+1}|_{2}=|\mathbf{w}^{\mathrm{new}}|_{2}$

16: Compute $q=\left(\frac{\mathrm{err}}{\mathrm{R_{TOL}}}\right)^{k_{I}}$,
$\mathrm{\mathrm{ratio}}=\frac{\Delta t^{\max}}{\Delta t}$,
$\mathrm{ratio}=\min(\mathrm{ratiomax},\max(\mathrm{ratiomin},\mathrm{ratio}))$

17: if $|\mathrm{ratio}-1|>\mathrm{ratiomin}$ then

18: $\Delta t=\mathrm{ratio}\cdot\Delta t$,
$\mathbf{w}_{z}^{n+1}=\mathrm{ratio}\cdot\mathbf{w}_{z}^{n+1}$

19: needNewMJ $\leftarrow$ true

20: end if

21:end while

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 $\mathbf{w}_{z}^{n+\gamma}$ can be simply taken as
$\mathbf{w}_{z}^{n}$, 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 $\mathbf{w}_{z,k+1}^{new}$ and
$\mathbf{w}_{k+1}^{new}$, a second Newton step is necessary within the
implicit stage (S2), the particularity of which is the evaluation of the
initial guess $\mathbf{w}_{z,0}^{new}$, by extrapolation of the derivatives of
the cubic Hermite interpolant from $t^{n}$ to $t^{n+\gamma}$ [29].

Algorithm 2 Adaptive Runge-Kutta time stepping: procedure TR-
BDF2($\mathbf{w}^{n}$,$\mathbf{w}_{z}^{n}$,$\Delta
t$,$\boldsymbol{\beta}_{\Delta t}$)

1:compute $\mathrm{sc}_{i}=\max(|\mathbf{w}^{n}|_{2},\eta)$

2:assemble $\mathbf{M}$ and $\mathbf{J}$ using
MJcontrol($\mathbf{w}^{n}$,$\mathbf{w}_{z}^{n}$,$\Delta
t$,$\boldsymbol{\beta}_{\Delta t}$) $\triangleright$ See Algorithm 4

3:set $t_{n+\gamma}=t+2\gamma\Delta t$,
$\mathbf{w}_{0}^{n+\gamma}=\mathbf{w}^{n}+2\gamma\mathbf{w}_{z}^{n}$,
$\mathbf{w}_{z,0}^{n+\gamma}=\mathbf{w}_{z}^{n}$ $\triangleright$ Stage 2

4:rejectS1 $\leftarrow$ false

5:set
$\mathrm{res}_{\mathrm{S_{2},min}}=100\epsilon\vvvert\mathbf{w}_{0}^{n+\gamma}\vvvert$

6:for $k=0,\ldots,K$ do

7: for $*\in\\{D,E\\}$ do

8: $\mathbf{g}_{k}^{*,n+\gamma}=-\frac{1}{\Delta
t}\mathbfsf{M}^{*}\mathbf{w}_{z,k}^{*,n+\gamma}+\left[-\mathbfsf{D}^{*}\mathbf{w}_{k}^{*,n+\gamma}+\mathbf{G}^{*}(\mathbf{w}_{k}^{*,n+\gamma})-\mathbfsf{C}_{u}^{*,n}\mathbf{w}_{k}^{*,n+\gamma}\right]$

9: if $\neg$ mjcontrol then

10: Update $\mathbfsf{M}^{*}=\mathbfsf{M}_{k}^{*,n+\gamma}$ and
$\mathbf{J}^{*}=\mathbf{J}_{k}^{*,n+\gamma}$

11: end if

12: solve $\delta_{k+1}^{n+\gamma}\mathbf{w}_{z}^{*}=(\frac{1}{\Delta
t}\mathbfsf{M}^{*}-\gamma\mathbfsf{J}^{*})^{-1}\mathbf{g}_{k}^{*,n+\gamma}$

13: compute
$\mathbf{w}_{z,k+1}^{*,n+\gamma}=\mathbf{w}_{z,k}^{*,n+\gamma}+\delta_{k+1}^{n+\gamma}\mathbf{w}_{z}^{*}$,
$\mathbf{w}_{k+1}^{*,n+\gamma}=\mathbf{w}_{k}^{*,n+\gamma}+\gamma\delta_{k+1}^{n+\gamma}\mathbf{w}_{z}^{*}$

14: update Robin transmission conditions $\mathcal{R}_{h}^{*}$ from (2.7)

15: end for

16:
NewtonConv($\mathbf{w}_{k+1}^{n+\gamma}$,$\delta_{k+1}^{n+\gamma}\mathbf{w}_{z}$,$\mathrm{res}_{\mathrm{S_{2},min}}$,$\boldsymbol{\beta}_{\Delta
t}$) $\triangleright$ See Alg. 3

17:end for

18:rejectS2 $\leftarrow$ false

19:if $\neg$ rejectS1 then$\triangleright$ Stage 3

20: compute
$\mathrm{sc}_{i}=\max(|\mathbf{w}^{n+\gamma}|_{2},\mathrm{sc}_{i})$

21:
$\mathbf{w}_{z,0}^{\mathrm{new}}=p_{31}\mathbf{w}_{z}^{n}+p_{32}\mathbf{w}_{z}^{n+\gamma}+p_{33}\mathbf{w}^{n+\gamma}-p_{33}\mathbf{w}^{n}$

22:
$\mathbf{w}_{0}^{\mathrm{new}}=\mathbf{w}^{n}+(1-b_{2}-\gamma)\mathbf{w}_{z}^{n}+b_{2}\mathbf{w}_{z}^{n+\gamma}+\gamma\mathbf{w}_{z,0}^{\mathrm{new}}$

23: redo steps 6-15 to retrieve $\mathbf{w}_{z,k+1}^{\mathrm{new}}$ and
$\mathbf{w}_{k+1}^{\mathrm{new}}$

24:end if

25:if rejectS1 $\lor$ rejectS2 then

26: if mjcontrol then

27: if Jcurrent $\land$ Mcurrent then

28: set $\Delta t^{\mathrm{old}}=\Delta t$, $\Delta t=\max(fac_{S1S2}\Delta
t^{\mathrm{old}},\Delta t^{\min})$, $\mathbf{w}_{z}^{n}=\frac{\Delta t}{\Delta
t^{\mathrm{old}}}\mathbf{w}_{z}^{n}$

29: needNewMJ $\leftarrow$ true

30: else

31: needNewJ $\leftarrow$ $\neg$Jcurrent, needNewM $\leftarrow$ $\neg$Mcurrent

32: end if

33: else

34: set $\Delta t^{\mathrm{old}}=\Delta t$, $\Delta t=\max(fac_{S1S2}\Delta
t^{\mathrm{old}},\Delta t^{\min})$, $\mathbf{w}_{z}^{n}=\frac{\Delta t}{\Delta
t^{\mathrm{old}}}\mathbf{w}_{z}^{n}$

35: end if

36:else

37: compute
$\mathrm{sc}_{i}=\max(|\mathbf{w}^{\mathrm{new}}|_{2},\mathrm{sc}_{i})$

38:
$\boldsymbol{e}_{w}^{1}=(\hat{b}_{1}-b_{1})\mathbf{w}_{z}^{n}+(\hat{b}_{2}-b_{2})\mathbf{w}_{z}^{n+\gamma}+(\hat{b}_{3}-b_{3})\mathbf{w}_{z}^{n+1}$

39: if $\neg$ mjcontrol then

40: Update $\mathbf{M}=\mathbf{M}^{\mathrm{new}}$ and
$\mathbf{J}=\mathbf{J}^{\mathrm{new}}$

41: end if

42: compute $\boldsymbol{e}_{w}^{2}=(\frac{1}{\Delta
t}\mathbfsf{M}-\gamma\mathbfsf{J})^{-1}\boldsymbol{e}_{w}^{1}$ and
$\mathrm{err}=\max(\vvvert\boldsymbol{e}_{w}^{2}\vvvert,\frac{\vvvert\boldsymbol{e}_{w}^{1}\vvvert}{16})$

43: if $\mathrm{err}>\mathrm{R_{TOL}}$ then

44: $\Delta t^{\mathrm{old}}=\Delta t$, $\Delta t=\max\left(\Delta
t^{\mathrm{old}}\max\left(fac_{\min},fac\left(\frac{\mathrm{R_{TOL}}}{err}\right)^{k_{I}}),\Delta
t^{\min}\right)\right)$, $\mathbf{w}_{z}^{n}=\frac{\Delta t}{\Delta
t^{\mathrm{old}}}\mathbf{w}_{z}^{n}$

45: needNewMJ $\leftarrow$ true

46: end if

47:end if

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
$\boldsymbol{v}$ of length $M$, here we adopt the rescaled maximum norm

$\vvvert\boldsymbol{v}\vvvert=\max_{i=1,\ldots,M}(v_{i}\cdot sc^{-1}_{i}),$

where $sc_{i}:=\max(sc_{i},|w_{k+1,i}|)$ is the local maximum of the solution
vector at the $k$-th Newton step, and it is reinitialised as
$sc_{i}=\max(|\mathbf{w}^{n}|_{2},\eta)$, where
$\eta=\frac{\mathrm{A_{TOL}}}{\mathrm{R_{TOL}}}$. 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 $\delta_{k+1}\mathbf{w}_{z}^{*}$. 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
$\|\delta_{k+1}\mathbf{w}_{z}\|\leq\theta\|\delta_{k}\mathbf{w}_{z}\|$, with
$\theta<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 $\frac{\theta}{1-\theta}\mathrm{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 $\kappa_{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
$\kappa_{N}\mathrm{TOL}_{N}<e(\delta_{k+1}\mathbf{w}_{z},\widehat{sc})\theta^{K-k}$
(cf. Line 26). Point (c) represents a rough estimate of the iteration error
expected after $K-1$ iterations.

Algorithm 3 Newton stopping criterion: procedure
NewtonConv($\mathbf{w}_{k+1}$,$\delta_{k+1}\mathbf{w}_{z}$,$\mathrm{res}_{\mathrm{\min}}$,$\boldsymbol{\beta}_{\Delta
t}$)

1:compute
$\widehat{\mathrm{sc}}_{i}=\max(\mathrm{sc}_{i},|\mathbf{w}_{k+1,i}|)$ and
$\mathrm{res}_{\mathrm{k+1}}(\widehat{\mathrm{sc}}_{i})=\vvvert\delta_{k+1}\mathbf{w}_{z}\vvvert$

2:if $\mathrm{res}_{\mathrm{k+1}}\leq\mathrm{res}_{\mathrm{\min}}$ then

3: break

4:else if $k=0$ then

5: if needNewRate $\land$ Stage 2 then

6: if mjcontrol then

7: needNewRate $\leftarrow$ false

8: end if

9: $\theta=0$

10: else

11:
$e(\delta_{k+1}\mathbf{w}_{z},\widehat{sc})=\frac{\theta}{1-\theta}\mathrm{res}_{\mathrm{k+1}}$

12: if $e(\delta_{k+1}\mathbf{w}_{z},\widehat{sc})\leq
0.1\kappa_{N}\mathrm{TOL}_{N}$ then

13: break

14: end if

15: end if

16:else if $\mathrm{res}_{\mathrm{k+1}}>0.9\mathrm{res}_{\mathrm{k}}$
then$\triangleright$ Divergence

17: rejectSX $\leftarrow$ true $\triangleright$ X stands for the stage

18: break

19:else

20: set
$\theta=\max\left(0.9\theta,\frac{\mathrm{res}_{\mathrm{k+1}}}{\mathrm{res}_{\mathrm{k}}}\right)$
and
$e(\delta_{k+1}\mathbf{w}_{z},\widehat{sc})=\frac{\theta}{1-\theta}\mathrm{res}_{\mathrm{k+1}}$

21: if
$e(\delta_{k+1}\mathbf{w}_{z},\widehat{sc})\leq\kappa_{N}\mathrm{TOL}_{N}$
then$\triangleright$ Convergence

22: break

23: else if $k=K$ then $\triangleright$ Max. iterations attained

24: rejectSX $\leftarrow$ true

25: break

26: else if
$\kappa_{N}\mathrm{TOL}_{N}<e(\delta_{k+1}\mathbf{w}_{z},\widehat{sc})\theta^{K-k}$
then$\triangleright$ Divergence

27: rejectSX $\leftarrow$ true

28: break

29: end if

30:end if

Should a stage be rejected, one moves to Line 25 in Alg. 2 and reduces the
timestep $\Delta t$ by a factor $fac_{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
$\mathbf{e}_{w}^{1}$ 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 $\Delta t\leq\Delta t^{\min}$.

Algorithm 4 Control of mass-stiffness assembly: procedure
MJcontrol($\mathbf{w}^{n}$,$\mathbf{w}_{z}^{n}$,$\Delta
t$,$\boldsymbol{\beta}_{\Delta t}$)

1:if needNewM$\land$mjcontrol then

2: compute $\mathbfsf{M}^{*,n}$, with $*\in\\{D,E\\}$

3: Mcurrent$\leftarrow$true, needNewM$\leftarrow$false,
needNewMJ$\leftarrow$true

4:end if

5:if needNewJ$\land$mjcontrol then

6: compute
$\mathbfsf{J}^{*,n}=-\mathbfsf{D}^{*}+\mathbfsf{J}_{\mathbf{G}^{*}}^{n}-\mathbfsf{C}_{u}^{*,n}$
with $*\in\\{D,E\\}$

7: Jcurrent$\leftarrow$true, needNewJ$\leftarrow$false,
needNewMJ$\leftarrow$true

8:end if

9:if needNewMJ$\land$mjcontrol then

10: compute $\frac{1}{\Delta t}\mathbfsf{M}^{*,n}-\gamma\mathbfsf{J}^{*,n}$
with $*\in\\{D,E\\}$

11: needNewMJ$\leftarrow$false, needNewRate$\leftarrow$true

12: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 exploited.

mjcontrol |  Switch to reuse Newton iteration matrices over two implicit stages
---|---
needNewM |  Need to evaluate mass matrix $\mathbf{M}^{*}$
needNewJ |  Need to evaluate Jacobian matrix $\mathbf{J}^{*}$
needNewMJ |  Need to evaluate mass-Jacobian matrix $\frac{1}{\Delta t}\mathbfsf{M}^{*}-\gamma\mathbfsf{J}^{*}$
Mcurrent |  Appropriate mass matrix $\mathbf{M}$ for the current timestep
Jcurrent |  Appropriate Jacobian matrix $\mathbf{J}$ for the current timestep
Table 3.2: Set of Booleans controlling the evaluation and assembly of mass and
stiffness matrices.

##  4 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).

$c_{f}^{\mathrm{D}}=150$, $c_{g}^{\mathrm{D}}=1$,
$\rho_{0}^{\mathrm{D}}=\rho_{2}^{\mathrm{D}}=\rho_{3}^{\mathrm{D}}=\rho_{5}^{\mathrm{D}}=1$,
$\rho_{1}^{\mathrm{D}}=0$, $\rho_{4}^{\mathrm{D}}=0.35$,
$M_{11}^{\mathrm{D}}=1$, $M_{22}^{\mathrm{D}}=30$,
---
$c_{f}^{\mathrm{E}}=20$, $c_{g}^{\mathrm{E}}=2$, $\rho_{0}^{\mathrm{E}}=2$,
$\rho_{2}^{\mathrm{E}}=\rho_{3}^{\mathrm{E}}=2$, $\rho_{5}^{\mathrm{E}}=1$,
$\rho_{1}^{\mathrm{E}}=0$, $\rho_{4}^{\mathrm{E}}=0.15$,
$M_{11}^{\mathrm{E}}=2$, $M_{22}^{\mathrm{E}}=10$,
$E^{\mathrm{D}}=1000$, $\nu^{\mathrm{D}}=0.475$, $J^{\mathrm{D}}=1$,
$K^{\mathrm{D}}=1$e5, $E^{\mathrm{E}}=10$, $\nu^{\mathrm{D}}=0.33$,
$J^{\mathrm{E}}=1$, $K^{\mathrm{E}}=1$e5, $\alpha^{\mathrm{E}}=2.5$
Table 4.1: Example 1. Parameters for the mechanochemical model using (2.3) and
(2.5).

#### Example 1: Experimental accuracy against manufactured solutions.

Let us consider a rectangular domain $\Omega=(0,1)\times(0,1.4)$ divided into
the dermis $\Omega^{\mathrm{D}}=(0,1)^{2}$ and epidermis
$\Omega^{\mathrm{E}}=\Omega\setminus\Omega^{\mathrm{D}}$ subdomains. The
initial position of the interface is characterised by the segment
$x\in(0,1),y=1$, and the external surface is $x\in(0,1),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:

$\tilde{\boldsymbol{w}}=\begin{pmatrix}1-\cos(2\pi x)\sin(3\pi y)\\\
1+\frac{1}{2}\cos(2\pi x)\sin(3\pi
y)\end{pmatrix},\quad\tilde{\boldsymbol{u}}=\begin{pmatrix}x(1-x)\cos(\pi
x)\sin(2\pi y)\\\ \sin(\pi x)\cos(\pi y)y^{2}(1-y)\end{pmatrix},$ (4.1)

and set
$\boldsymbol{u}^{\mathrm{D}}=\tilde{\boldsymbol{u}}|_{\Omega^{\mathrm{D}}}$,
$\boldsymbol{u}^{\mathrm{E}}=\tilde{\boldsymbol{u}}|_{\Omega^{\mathrm{E}}}$,
$\boldsymbol{w}^{\mathrm{D}}=\tilde{\boldsymbol{w}}|_{\Omega^{\mathrm{D}}}$,
$\boldsymbol{w}^{\mathrm{E}}=\tilde{\boldsymbol{w}}|_{\Omega^{\mathrm{E}}}$,
and $p^{*}=-\lambda^{*}\operatorname*{div}\tilde{\boldsymbol{u}}^{*}$. The
functions in (4.1) are used to construct the Dirichlet datum for the
displacements on $\partial\Omega\setminus\Gamma$ as well as the flux
conditions imposed on the morphogens at the boundary $\partial\Omega$.
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.

Space convergence
---
DoF | $h$ | $e_{1}(\boldsymbol{w}^{\mathrm{D}})$ | $r_{1}(\boldsymbol{w}^{\mathrm{D}})$ | $e_{1}(\boldsymbol{w}^{\mathrm{E}})$ | $r_{1}(\boldsymbol{w}^{\mathrm{E}})$ | $e_{0}(p^{\mathrm{D}})$ | $r_{0}(p^{\mathrm{D}})$ | $e_{0}(p^{\mathrm{E}})$ | $r_{0}(p^{\mathrm{E}})$ |
48 | 0.5207 | 5.6275 | – | 3.5698 | – | 1175.8 | – | 2.2114 | – |
108 | 0.2828 | 4.1858 | 0.5794 | 2.6495 | 0.4886 | 444.12 | 1.9061 | 1.2433 | 0.9437 |
280 | 0.1736 | 2.6223 | 0.7956 | 1.8701 | 0.7133 | 129.18 | 2.1008 | 0.6340 | 1.3790 |
900 | 0.0889 | 1.4380 | 0.9447 | 0.9974 | 0.9398 | 37.750 | 1.9344 | 0.2178 | 1.5973 |
3264 | 0.0432 | 0.7485 | 0.9844 | 0.4852 | 0.9978 | 11.335 | 1.8138 | 0.0801 | 1.3849 |
12276 | 0.0218 | 0.3829 | 0.9890 | 0.2480 | 0.9788 | 3.6282 | 1.6805 | 0.0409 | 0.9624 |
47320 | 0.0110 | 0.1970 | 0.9695 | 0.1326 | 0.9217 | 1.4185 | 1.3701 | 0.0252 | 0.9748 |
186276 | 0.0055 | 0.1065 | 0.9928 | 0.0684 | 0.9605 | 0.7001 | 0.9599 | 0.0141 | 0.9635 |
DoF | $h$ | $e_{0}(\boldsymbol{u}^{\mathrm{D}})$ | $r_{0}(\boldsymbol{u}^{\mathrm{D}})$ | $e_{1}(\boldsymbol{u}^{\mathrm{D}})$ | $r_{1}(\boldsymbol{u}^{\mathrm{D}})$ | $e_{0}(\boldsymbol{u}^{\mathrm{E}})$ | $r_{0}(\boldsymbol{u}^{\mathrm{E}})$ | $e_{1}(\boldsymbol{u}^{\mathrm{E}})$ | $r_{1}(\boldsymbol{u}^{\mathrm{E}})$ | Iter
48 | 0.5207 | 0.0605 | – | 1.2176 | – | 0.0423 | – | 0.5038 | – | 4
108 | 0.2828 | 0.0255 | 1.6925 | 0.6149 | 1.3372 | 0.0192 | 1.2900 | 0.3408 | 0.6409 | 4
280 | 0.1736 | 0.0081 | 1.9584 | 0.2545 | 1.5009 | 0.0093 | 1.4921 | 0.2226 | 0.8723 | 4
900 | 0.0889 | 0.0024 | 1.9235 | 0.1201 | 1.1808 | 0.0029 | 1.7335 | 0.1131 | 1.0121 | 4
3264 | 0.0432 | 0.0006 | 1.9717 | 0.0596 | 1.0563 | 0.0011 | 1.2942 | 0.0537 | 1.0313 | 4
12276 | 0.0218 | 0.0002 | 1.9373 | 0.0299 | 1.0180 | 0.0006 | 1.9020 | 0.0275 | 0.9765 | 4
47320 | 0.0110 | 3.61e-5 | 1.9734 | 0.0150 | 1.0061 | 3.70e-5 | 1.9258 | 0.0151 | 0.9848 | 4
186276 | 0.0055 | 9.02e-6 | 1.8221 | 0.0075 | 1.0020 | 9.12e-6 | 1.8300 | 0.0065 | 0.9730 | 4

Time convergence
---
$\Delta t$ | $e_{\Delta t}(\boldsymbol{w}^{\mathrm{D}})$ | $r_{\Delta t}(\boldsymbol{w}^{\mathrm{D}})$ | $e_{\Delta t}(\boldsymbol{w}^{\mathrm{E}})$ | $r_{\Delta t}(\boldsymbol{w}^{\mathrm{E}})$ | avg(iter)
0.0100 | 0.2002 | – | 0.0881 | – | 3
0.0050 | 0.0982 | 1.1294 | 0.0292 | 1.1752 | 3
0.0025 | 0.0230 | 1.9301 | 0.0094 | 1.7388 | 3
0.0013 | 0.0056 | 1.9614 | 0.0024 | 1.7929 | 3
0.0006 | 0.0016 | 1.9811 | 0.0006 | 1.9952 | 3
0.0002 | 0.0004 | 1.9937 | 0.0002 | 1.9874 | 3

Table 4.2: Example 1. Spatial and temporal error history associated to the
discretisation of the model problem. Errors on the bottom sub-table were
obtained on a fine mesh of size $h=0.0013$, and computed until a final time
$T=0.02$.

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}(s^{*}):=\|s^{*}-s^{*}_{h}\|_{0,\Omega^{*}},\quad
e_{1}(s^{*}):=\|s^{*}-s^{*}_{h}\|_{1,\Omega^{*}},\quad
r_{k}(s^{*})=\frac{\log(e_{k}(s^{*})/\widehat{e_{k}}(s^{*}))}{\log(h/\widehat{h})},\quad
k=0,1,$

where $e_{k}$ and $\widehat{e_{k}}$ denote errors computed on two consecutive
meshes of sizes $h$ and $\widehat{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 $\mathrm{TOL}_{N}=1$e-5, and the stopping criterion is
based on the $\ell^{\infty}-$norm 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

$\boldsymbol{w}=\tilde{\boldsymbol{w}}\exp(-kt),$

with $k=\log(1/2)$. Errors and convergence rates are measured and denoted as
follows

$e_{\Delta
t}(s^{*}):=\sum_{n=1}^{N}\|s^{*}(t^{n})-s^{*}_{h}(t^{n})\|_{0,\Omega^{*}},\quad
r_{\Delta t}(s^{*})=\frac{\log(e_{\Delta t}(s^{*})/\widehat{e_{\Delta
t}}(s^{*}))}{\log(\Delta t/\widehat{\Delta t})},$

where $e_{\Delta t},\widehat{e_{\Delta t}}$ denote errors computed using two
consecutive timesteps $\Delta t$ and $\widehat{\Delta 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]).

$\Delta t^{max}$ | 2000 | $\mathrm{R_{TOL}}$ | $10^{-6}$ | $\mathrm{A_{TOL}}$ | $10^{-3}$
---|---|---|---|---|---
$\mathrm{TOL}_{N}$ | $10^{-6}$ | $\kappa_{N}$ | 0.5 | $K$ | 10
$fac_{S1S2}$ | 0.3 | $fac_{\min}$ | 0.1 | $fac$ | $\sqrt[3]{0.25}$
$\mathrm{ratiomin}$ | 0.2 | $\mathrm{ratiomax}$ | 5.0 | $\mathrm{ratio}$ | $\sqrt[3]{0.25}$
$k_{I}$ | $\frac{1}{3}$ | $\epsilon$ | $10^{-10}$ |  |
Table 4.3: Example 2. Adaptive timestep controller parameters.

#### 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
$\mathbf{M}^{*}_{ij}\neq 0$ $\forall 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 $\rho_{0}=0$, $\rho_{1}=\rho_{2}=\rho_{4}=\rho_{5}=1$,
$\rho_{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^{\mathrm{D}}=K^{\mathrm{E}}=1$. The domains of interest are the blocks
$\Omega^{\mathrm{D}}=(0,50)^{2}$ and
$\Omega^{\mathrm{E}}=(0,50)\times(50,75)$, which we discretise into
unstructured meshes with 16360 and 8032 triangular elements, respectively.
Initial concentrations are set according to

$\tilde{w}_{1}=\frac{\rho_{2}}{\rho_{3}}\biggl{(}\rho_{0}+\rho_{1}\frac{\rho_{5}}{\rho_{4}}\biggr{)},\quad\tilde{w}_{2}=\frac{\rho_{4}}{\rho_{5}}\tilde{w}_{1}^{2},\quad\boldsymbol{w}^{*}_{0}(\boldsymbol{x})=\begin{pmatrix}\tilde{w}_{1}(1+\eta(\boldsymbol{x}))\\\
\tilde{w}_{2}\end{pmatrix},$

where $\eta(\boldsymbol{x})$ 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
Figure 4.1.

Figure 4.1: Example 2. Species concentrations $w_{1}^{*}$, $w_{2}^{*}$
computed at $T=2000$.

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 $\Delta 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
$\texttt{mjcontrol}=\textbf{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.

Figure 4.2: Example 2. Time step evolution and number of Newton iterations
needed to converge on stages S1 and S2, and setting mjcontrol on and off.

In addition, we test our algorithm for both linear and nonlinear cross-
diffusion terms defined as

$\mathbf{M}_{lin}^{*}(\boldsymbol{w}^{*})=\begin{pmatrix}M^{*}_{11}&M^{*}_{12}\\\
M^{*}_{21}&M^{*}_{22}\end{pmatrix}$

or

$\mathbf{M}_{\mathrm{nl}}^{*}(\boldsymbol{w}^{*})=\begin{pmatrix}M^{*}_{11}(1+\eta_{11}^{*,w_{1}}w_{1}^{*}+\eta_{11}^{*,w_{2}}w_{2}^{*})&0\\\
0&M^{*}_{22}(1+\eta_{22}^{*,w_{1}}w_{1}^{*}+\eta_{22}^{*,w_{2}}w_{2}^{*})\end{pmatrix},$

where $M_{ij}^{*}$, $\eta_{ij}^{*,w_{k}}$ ($i,j,k=\\{1,2\\}$) are real
constants. For these simulations we impose $M^{*}_{11}=M^{*}_{1}$,
$M^{*}_{22}=M^{*}_{2}$. 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 cross-diffusion models, we vary spatially the values $[M_{12}^{*}$,
$M_{21}^{*}]\in[0,1]\times[0,15]$, and
$\eta_{11}^{*,w_{1}}=\eta_{22}^{*,w_{1}}=\eta_{1}$ and
$\eta_{11}^{*,w_{2}}=\eta_{22}^{*,w_{2}}=\eta_{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 ).

Figure 4.3: Example 2, 2D case. Concentration of species $w_{1}^{*}$. Linear
cross-diffusion (top row) using $(M_{12}^{*},M_{21}^{*})=(0.5,0)$, $(0,0.5)$
and $(0.35,0.05)$ (left, centre, right, respectively). The bottom panels show
the case of nonlinear cross-diffusion using
$(\eta_{11}^{*,w_{1}},\eta_{11}^{*,w_{2}},\eta_{22}^{*,w_{1}},\eta_{22}^{*,w_{2}}=(10^{-1},10^{-1},10^{-1},10^{-1})$,
$(1,0,0,1)$ and $(1,0,1,0)$ (left, centre, right, respectively).

Figure 4.4: Example 2. Space parameter plots produced with linear (left) and
nonlinear (right) cross-diffusion terms.

Figure 4.5: Example 2, 3D case. Concentration of species $w_{1}^{*}$. Linear
cross-diffusion (top row) using $(M_{12}^{*},M_{21}^{*})=(0.5,0)$, $(0,0.5)$
and $(0.35,0.05)$ (left, centre, right, respectively). The bottom panels show
the case of nonlinear cross-diffusion using
$(\eta_{11}^{*,w_{1}},\eta_{11}^{*,w_{2}},\eta_{22}^{*,w_{1}},\eta_{22}^{*,w_{2}})=(10^{-1},10^{-1},10^{-1},10^{-1})$,
$(1,0,0,1)$ and $(1,0,1,0)$ (left, centre, right, respectively).

Figure 4.6: Example 3. Morphogen concentrations $w_{1}^{*}$, $w_{2}^{*}$ (left
and centre) and displacement magnitude (right), computed at $T=2000$ and
plotted on the deformed configuration.

#### 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^{\mathrm{D}}=1000$, $E^{\mathrm{E}}=250$, $\nu^{\mathrm{D}}=0.475$,
$\nu^{\mathrm{E}}=0.3$. Displacements are prescribed on the boundary
$\partial\Omega\setminus\Gamma$ according to
$\boldsymbol{u}=(\frac{1}{2}\cos(5\pi(x-y)/2),\frac{3}{4}\sin(5\pi(x+y)/2))^{T}$.
Additionally, we impose $c_{f}^{\mathrm{D}}=150$, $c_{f}^{\mathrm{E}}=20$ and
$c_{g}^{*}=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 $\mathrm{A_{TOL}}=10^{-6}$ and
$\mathrm{TOL}_{N}=\mathrm{R_{TOL}}=10^{-4}$. The process is simulated until
$T=2000$.

Figure 4.6 presents the final patterns generated by the morphogen
concentrations and also the zones of high displacement. On the bottom layer
the pattern is similar to the one in 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 $\Omega^{\mathrm{D}}$.

Figure 4.7: Example 3. Concentration of species $w_{1}^{*}$ and displacement
magnitude $|u^{*}|$. Influence of ADR parameters ($1^{\text{st}}$ &
$2^{\text{nd}}$ rows) using $M_{1}^{*}=0.5$, $M_{1}^{*}=1.5$,
$\rho_{2}^{*}=0.5$ and $\rho_{3}^{*}=0.45$ (from left to right respectively);
and of coupling parameters ($3^{\text{rd}}$ & $4^{\text{th}}$ rows) using
$c_{f}^{*}=150$, $c_{f}^{*}=20$, $c_{g}^{*}=0.5$, $c_{g}^{*}=1.5$ (from left
to right respectively).

We next test the influence of the reaction and mechanical coupling parameters
$(c_{f}^{*},c_{g}^{*})$ on the final pattern, and typical results are depicted
in Figure 4.7. Perturbations of the ADR coefficients 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 $\Sigma$. 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}_{Newt}$ represents the number of iterations
corresponding to _accepted_ timesteps. For this case we observe that $\Delta
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 $\nu^{*}=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 converge in less than 4
iterations in both stages of the RK method.

Figure 4.8: Example 3. Top: Timestep evolution and Newton iteration count at
stages S1 and S2 for different mechanical properties and acceleration
constants. Species concentration $w_{1}^{*}$ using $\nu^{*}=0.3$,
$\nu^{*}=0.475$, $J^{D}=10$ and $J^{D}=0.1$ (bottom, from left to right,
respectively).

#### 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
$\Omega^{\mathrm{D}}=(0,50)^{3}$ and
$\Omega^{\mathrm{E}}=(0,50)^{2}\times(50,75)$, 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 $\mathrm{A_{TOL}}=10^{-6}$
and $\mathrm{TOL}_{N}=\mathrm{R_{TOL}}=10^{-3}$. The process is run until
$T=300$.

$c_{f}^{\mathrm{D}}=20$, $c_{g}^{\mathrm{D}}=1$, $r_{1}=0.25$, $r_{0}=1$,
$r_{2}=20$, $r_{3}=5$, $r_{4}=5$, $M_{11}^{\mathrm{D}}=M_{22}^{\mathrm{D}}=1$,
$M_{33}^{\mathrm{D}}=35$, $M_{44}^{\mathrm{D}}=10$,
---
$c_{f}^{\mathrm{E}}=20$, $c_{g}^{\mathrm{E}}=1$ $r_{5}=5.0$, $r_{6}=20.0$,
$r_{7}=5.0$, $M_{11}^{\mathrm{E}}=M_{22}^{\mathrm{E}}=1$,
$M_{33}^{\mathrm{E}}=35$, $M_{44}^{\mathrm{E}}=10$,
$E^{*}=1000$, $\nu^{*}=0.475$, $J^{*}=K^{*}=1$, $\alpha^{\mathrm{E}}=9.5$
Table 4.4: Example 4. Parameters for the mechanochemical model using (2.4) and
(2.5).

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 $w_{2}^{*}$ (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 $\Omega$, 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 $w_{1}^{*}$, 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 $w_{3}^{*}$ and $w_{4}^{*}$
we observe a similar behaviour, but having a larger variation on the $z$-axis
for a specific layer ($\Omega^{\mathrm{D}}$ for $w_{4}$, $\Omega^{\mathrm{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 degradation, as one moves away from that layer.

Figure 4.9: Example 4, 3D case. Concentration of species $w_{1}^{*}$,
$w_{4}^{*}$ (left and and centre) and displacement magnitude (right), computed
a $T=300$ and plotted on the deformed configuration (scaled 10 times for
better visualisation). For the displacement magnitude, we cut the domain in
half to see the internal displacement.

##  5 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 advection-diffusion-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.

## References

  * [1] L. Ahtiainen, S. Lefebvre, P.H. Lindfors, E. Renvoisé, V. Shirokova, M.K. Vartiainen, I. Thesleff, and M.L. Mikkola, Directional cell migration, but not proliferation, drives hair placode morphogenesis. Develop. Cell (2014), 28(5):588–602.
  * [2] R. Allena, J.J. Muñoz, and D. Aubry, Diffusion-reaction model for Drosophila embryo development. Comput. Methods Biomech. Biomed. Engrg. (2013), 16(3):235–248.
  * [3] A. Aman and T. Piotrowski, Cell migration during morphogenesis, Dev. Biol. (2010), 341(1):20–33.
  * [4] B. Andreianov, M. Bendahmane, A. Quarteroni, and R. Ruiz-Baier, Solvability analysis and numerical approximation of linearized cardiac electromechanics, Math. Models Meth. Appl. Sci. (2015), 25(05):959–993.
  * [5] D.N. Arnold, F. Brezzi, and M. Fortin, A stable finite element for the Stokes equations. Calcolo (1984), 21:337–344.
  * [6] A. Badugu, C. Kraemer, P. Germann, D. Manshykau, and D. Iber, Digit patterning during limb development as a result of the BMP-receptor interaction, Sci Rep. (2012), 2:991.
  * [7] R.E. Bank, W.M. Coughran Jr, W. Fichtner, E.H. Grosse, D.J. Rose, and R.K. Smith, Transient simulation of silicon devices and circuits, IEEE Trans. Comput.-Aided Design (1985), 4:436–451.
  * [8] R. Becker, E. Burman, and P. Hansbo, A Nitsche extended finite element method for incompressible elasticity with discontinuous modulus of elasticity, Comput. Methods Appl. Mech. Engrg. (2009), 198:3352–3360.
  * [9] F. BenBelgacem, L.K. Chilton, and P. Seshaiyer, The $hp$-mortar finite-element method for the mixed elasticity and Stokes problems, Comput. Math. Appl. (2003), 46:35–55.
  * [10] L. Bonaventura and A. Della Rocca, Unconditionally strong stability preserving extensions of the TR-BDF2 method, J. Sci. Comput. (2017), 70:859–895.
  * [11] M. Buck, O. Iliev, and H. Andrä, Multiscale finite element coarse spaces for the application to linear elasticity, Central Eur. J. Math. (2013), 11(4):680–701.
  * [12] W. Chen, X. Xu, and S. Zhang, On the optimal convergence rate of a Robin-Robin domain decomposition method, J. Comput. Math. (2014), 32(4):456–475.
  * [13] A. Chertock, K. Fellner, A. Kurganov, A. Lorz, and P.A. Markowich, Sinking, merging and stationary plumes in a coupled chemotaxis-fluid model: a high-resolution numerical approach, J. Fluid Mech. (2012), 694:155–190.
  * [14] C. Cherubini, S. Filippi, A. Gizzi, and R. Ruiz-Baier, A note on stress-driven anisotropic diffusion and its role in active deformable media, J. Theor. Biol. (2017), 430(7):221–228.
  * [15] G.C. Cruywagen and J.D. Murray, On a tissue interaction model for skin pattern formation, J. Nonl. Sci. (1992), 2(2):217–240.
  * [16] R. Dillon and H.G. Othmer, A mathematical model for outgrowth and spatial patterning of the vertebrate limb bud, J. Theor. Biol. (1999), 197:295–330.
  * [17] N. Di-Poï and M.C. Milinkovitch, The anatomical placode in reptile scale morphogenesis indicates shared ancestry among skin appendages in amniotes, Sci. Adv. (2016), 2(6).
  * [18] V. Dolean, P. Jolivet, and F. Nataf, An Introduction to Domain Decomposition Methods. Algorithms, Theory, and Parallel Implementation, SIAM Philadelphia, PA, USA (2015).
  * [19] J.D. Edwards, J.E. Morel, and D.A. Knoll, Nonlinear variants of the TR/BDF2 method for thermal radiative diffusion, J. Comput. Phys. (2011), 230:1198–1214.
  * [20] F. Hecht, New development in FreeFem++, J. Numer. Math. (2012), 20(3-4):251–265.
  * [21] M.J. Gander, F. Magoules, and F. Nataf, Optimized Schwarz methods without overlap for the Helmholtz equation, SIAM J. Sci. Comput. (2002), 24:38–60.
  * [22] G.N. Gatica, B. Gomez-Vargas, and R. Ruiz-Baier, Analysis and mixed-primal finite element discretisations for stress-assisted diffusion problems, Comput. Methods Appl. Mech. Engrg. (2018), 337.
  * [23] G.N. Gatica and W.L. Wendland, Coupling of mixed finite elements and boundary elements for a hyperelastic interface problem, SIAM J. Numer. Anal. (1997), 34(6):2335–2356.
  * [24] A. Gierer and H. Meinhardt, A theory of biological pattern formation, Kybernetik (1972), 12:30–39.
  * [25] J.D. Glover, K.L. Wells, F. Matthäus, K.J. Painter, W. Ho, J. Riddell, J.A. Johansson, M.J. Ford, C.A.B. Jahoda, V. Klika, R.L. Mort, and D.J. Headon, Hierarchical patterning modes orchestrate hair follicle morphogenesis, PLOS Bio. (2007), 15(7):1–31.
  * [26] P. Gosselet and C. Rey, Non-overlapping domain decomposition methods in structural mechanics, Arch. Comput. Methods Engrg. (2007), 13(4):515–572.
  * [27] E. Hairer, S.P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I, Nonstiff Problems.  Springer Series in Computational Mathematics. Springer New York, 2008.
  * [28] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems.  Springer Series in Computational Mathematics. Springer New York, 2010.
  * [29] M.E. Hosea and L.F. Shampine, Analysis and implementation of TR-BDF2, Appl. Num. Math. (1996), 20(1):21–37.
  * [30] T. Lecuit and P.-F. Lenne, Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis, Nat. Rev. Mol. Cell. Biol. (2007), 8(8):633–644.
  * [31] P. Lenarda, M. Paggi, and R. Ruiz Baier, Partitioned coupling of advection-diffusion-reaction systems and Brinkman flows, J. Comput. Phys. (2017), 344:281–302.
  * [32] P.L. Lions, On the Schwarz alternating method. III: a variant for non-overlapping subdomains, In: T.F. Chan, R. Glowinski, J. Périaux, and O. Widlund, editors, Third International Symposium on Domain Decomposition Methods for Partial Differential Equations, held in Houston, Texas, March 20-22, 1989, Philadelphia, PA, 1990. SIAM.
  * [33] R.T. Liu, S.S. Liaw, and P.K. Maini, Two-stage Turing model for generating pigment patterns on the leopard and the jaguar, Phys. Rev. E (2006), 74(1):011914.
  * [34] D. Manoussaki, A mechanochemical model of angiogenesis and vasculogenesis, ESAIM: Math. Model. Numer. Anal. (2003), 37(04):581–599.
  * [35] M. Mercker, F. Brinkmann, A. Marciniak-Czochra, and T. Richter, Beyond Turing: mechanochemical pattern formation in biological tissues, Biol. Dir. (2016), 11(1):22–37.
  * [36] S.E. Millar, Molecular mechanisms regulating hair follicle development, J. Inv. Derm. (2002), 118(2):216–225.
  * [37] S.A. Montandon, A.C. Tzika, A.F. Martins, B. Chopard, and M.C. Milinkovitch, Two waves of anisotropic growth generate enlarged follicles in the spiny mouse, EvoDevo (2014), 5:33.
  * [38] P. Moreo, E.A. Gaffney, J.M. García-Aznar, and M. Doblaré, On the modelling of biological patterns with mechanochemical models: Insights from analysis and computation, Bull. Math. Biol. (2010), 72:400–431.
  * [39] J.D. Murray and G.F. Oster, Generation of biological pattern and form, IMA J. Math. Appl. Med. Biol. (1984), 1(1):51–75.
  * [40] J.D. Murray, P.K. Maini, and R.T. Tranquillo, Mechanochemical models for generating biological pattern and form in development, Phys. Reports (1988), 171(2):59–84.
  * [41] J.D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications.  Interdisciplinary Applied Mathematics. Springer New York, 2006.
  * [42] T. Oliver, J. Lee, and K. Jacobson, Forces exerted by locomoting cells, Cell Biol. (1994), 5(3):139–147.
  * [43] G.F. Oster, J.D. Murray, and A.K. Harris, Mechanical aspects of mesenchymal morphogenesis, Development (1983), 78(1):83–125.
  * [44] K.J. Painter, G.S. Hunt, K.L. Wells, J.A. Johansson, and D.J. Headon, Towards an integrated experimental–theoretical approach for assessing the mechanistic basis of hair and feather morphogenesis, Interface Focus (2012), 2(4):433–450.
  * [45] M. Pawlaczyk, M. Lelonkiewicz, and M. Wieczorowski, Age-dependent biomechanical properties of the skin, Postepy. Dermatol. Alergol. (2013), 30(5):302–306.
  * [46] M. Ptashnyk and B. Seguin, Homogenization of a system of elastic and reaction-diffusion equations modelling plant cell wall biomechanics, ESAIM: Math. Model. Numer. Anal. (2016), 50(2):593–631.
  * [47] A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential Equations, Oxford Science Publications, Oxford (1999).
  * [48] C. Remani, Numerical Methods for Solving Systems of Nonlinear Equations, 2012, Lakehead University, Thunder Bay, Ontario, Canada.
  * [49] R. Ruiz-Baier, Primal-mixed formulations for reaction-diffusion systems on deforming domains, J. Comput. Phys. (2015), 299:320–338.
  * [50] S. Sick, S. Reinker, J. Timmer, and T. Schlake, WNT and DKK Determine Hair Follicle Spacing Through a Reaction-Diffusion Mechanism, AAAS (2006), 314(5804):1447–1450.
  * [51] L.J. Shaw and J.D. Murray, Analysis of a model for complex skin patterns, SIAM J. App. Math. (1990), 50(2):628–648.
  * [52] H. Shoji, A. Mochizuki, Y. Iwasa, M. Hirata, T. Watanabe, S. Hioki, and S. Kondo, Origin of directionality in the fish stripe pattern, Dev. Dyna. (2003), 226(4):627-633.
  * [53] A.E. Shyer, A.R. Rodrigues, G.G. Schroeder, E. Kassianidou, S. Kumar, and R.M. Harland, Emergent cellular self-organization and mechanosensation initiate follicle pattern in the avian skin, Science (2017), 357(6353):811–815.
  * [54] G. Søderlind, Automatic control and adaptive time-stepping, Num. Alg. (2002), 31(1):281–310.
  * [55] L. Qin and X. Xu, Optimized Schwarz methods with Robin transmission conditions for parabolic problems, SIAM J. Sci. Comput. (2008), 31:608–623.
  * [56] A. Quarteroni, T. Lassila, S. Rossi, and R. Ruiz-Baier, Integrated Heart – Coupling multiscale and multiphysics models for the simulation of the cardiac function, Comput. Methods Appl. Mech. Engrg. (2017), 314:345–407.
  * [57] A. Tosin, Mechanics and chemotaxis in the morphogenesis of vascular networks, B. Math. Biol. (2006), 68(7):1819–1836.
  * [58] L. Tranqui and P. Tracqui, Mechanical signalling and angiogenesis. The integration of cell–extracellular matrix couplings, C.R. Acad. Sci. (2000), 323(1):31-47.
  * [59] A.M. Turing, The chemical basis of morphogenesis, Phil. Trans. R. Soc. London B (1990), 237:37–72.
  * [60] B.L. Vaughan, R.E. Baker, D. Kay, and P.K. Maini, A modified Oster-Murray-Harris mechanical model of morphogenesis, SIAM J. Appl. Math. (2013), 73:2124–2142.
  * [61] M.A. Wyczalkowski, Z. Chen, B.A. Filas, V.D. Varner, and L.A. Taber, Computational models for mechanics of morphogenesis, Birth Defects Res. (2012), 92(6):132–152.
  * [62] B. Wang, K. Xia, and G.W. Wei, Matched interface and boundary method for elasticity interface problems, J. Comput. Appl. Math. (2015), 285:203–225.