Post-Turing tissue pattern formation: Advent of mechanochemistry

Chemical and mechanical pattern formation is fundamental during embryogenesis and tissue development. Yet, the underlying molecular and cellular mechanisms are still elusive in many cases. Most current theories assume that tissue development is driven by chemical processes: either as a sequence of chemical patterns each depending on the previous one, or by patterns spontaneously arising from specific chemical interactions (such as “Turing-patterns”). Within both theories, mechanical patterns are usually regarded as passive by-products of chemical pre-patters. However, several experiments question these theories, and an increasing number of studies shows that tissue mechanics can actively influence chemical patterns during development. In this study, we thus focus on the interplay between chemical and mechanical processes during tissue development. On one hand, based on recent experimental data, we develop new mechanochemical simulation models of evolving tissues, in which the full 3D representation of the tissue appears to be critical for obtaining a realistic mechanochemical behaviour. The presented modelling approach is flexible and numerically studied using state of the art finite element methods. Thus, it may serve as a basis to combine simulations with new experimental methods in tissue development. On the other hand, we apply the developed approach and demonstrate that even simple interactions between tissue mechanics and chemistry spontaneously lead to robust and complex mechanochemical patterns. Especially, we demonstrate that the main contradictions arising in the framework of purely chemical theories are naturally and automatically resolved using the mechanochemical patterning theory.


Structural mechanics and notation
In the following, structural dynamics are mainly expressed in the Lagrangian or particle-centered framework. Let X be a particle in the undeformed initial configuration Ω and x = x(X, t) be its current position in the deformed one Ω(t) at time t. Thus the deformation is defined by χ : Ω × I → Ω(t) (with x = χ(X, t)) and is assumed to be invertible as well as continuously differentiable in space and time. Furthermore, the vector field U(X, t) = x(X, t) − X in Lagrangian description, or u(x, t) = x − X(x, t) in the Eulerian one, joining these positions is called displacement. This directly yields the relation u(x, t) = U(χ −1 (x, t), t) = U(X, t) which allows us to state all results in terms of u in the following. The same holds for the velocity field V(X, t). Denoting the material line element to a material curve by dX and a spatial line element or spatial tangent vector to a spatial curve by dx, then we obtain the important relation for the local change of relative positions by dx = F(X, t)dX, where F is the deformation gradient which is defined as F(X, t) := ∇χ(X, t) = ∇u(x, t) + I, J(X, t) := det(F(X, t)), with the identity matrix I. J represents the relative change of the local tissue volume and F plays an important role when we consider active deformations such as growth (c.f., following section). Based on F, the following strain tensors can be derived: The right hand Cauchy-Green tensor as well as the Green-Lagrange strain tensor both representing a measure of change of volume respectively shape. For further details on strain tensors in the material configuration we refer to Holzapfel [1], section 2.5.
Furthermore, we denote an infinitesimal force acting on a surface element by df and surface elements by dS and ds and normal directions by N and n in the reference and the current configuration, respectively. Then, we claim that there exist unique second order tensor fields σ(x, t), P(X, t) (Cauchy's stress theorem) such that df = σ(x, t)nds = P(X, t)NdS, where σ(x, t) is the Cauchy stress tensor and P(X, t) the first Piola-Kirchhoff stress tensor. Finally, we introduce Σ(X) = F −1 (x, t)P(X, t) which is the transformation of P to the reference configuration and is called second Piola-Kirchhoff stress tensor and expresses the stress solely in the reference configuration. All three stress tensors are related via P = JσF −T = FΣ; for more details on the concept of stress we refer to textbooks, e.g. Holzapfel [1], chapter 3.

Structural model equations
The shape of structural equations is determined by the two principles of conservation of mass and momentum. These balance principles are naturally given in the current configuration Ω(t). In strong formulation, they read for current and initial mass distributions ρ(x, t) and ρ 0 (X), external forces f (x, t) and surface forces σ(x, t)n which can be expressed with the stress tensor σ as discussed before (Eq. (10)). The divergence comes up by the application of the Gauss's theorem to Cauchy's first equation of motion. As mentioned before, it is usually more convenient to express continuum mechanics for elastic solids in the reference configuration. Then, Cauchy's first equation of motion in Lagrangian coordinates holds pointwise in the reference configuration Ω and reads We want to emphasise that this equation was derived from the physical principles of conservation of mass and momentum only. Especially, they are valid for all materials at all times and involve no additional modelling assumption. If we want to predict the material response of particular materials such as biological tissue, we additionally need to specify a material law e.g. relating the stress to the deformation gradient F. In particular, we restrict ourselves to hyperelastic materials where we assume the existence of a scalar-valued stored-energy function Ψ = Ψ(F), see e.g. Ciarlet [2] or Holzapfel [1]. There are many equivalent forms of the stored-energy function, especially it can be expressed via Ψ(F) = Ψ(E) = Ψ(C). This implies, that the material is homogeneous, since Ψ depends on the deformation gradient alone and not on the current position X. Furthermore, it is assumed that our biological tissue is isotropic and isothermal. This means, that the material response is the same in all directions and for all temperatures. We obtain relations of the form Obviously, assuming isotropy is a strong idealization: Each biological cell has a complex network of proteins called cytoskeleton, which gives a cell its shape and mechanical resistence to deformation. On the other hand, assuming isothermal behaviour is a reasonable assumption since temperature does not change significantly during biological processes such as pattern formation.
In the following, we use the simple nonlinear Saint-Venant-Kirchhoff model for compressible, hyperelastic materials to model the biological tissue. It states a linear relation between stress and strain, called Hooke's law. The strain energy function for the Saint-Venant-Kirchhoff model is given by which provides a relation between the second Piola-Kirchhoff stress tensor Σ and the Green-Lagrangian strain E tensor, namely where µ and λ are the Lamé constants. Values of the Lamé constants are usually given in terms of the Young's modulus E which is a measure of the stiffness (relation of stress to strain along an axis) of a solid material and Poisson's ratio ν which is the signed ratio of transverse to axial strain. Both measures can be obtained by the conversion formulas Note, that the Poisson's ratio for typical isotropic, compressible materials has a ratio of 0.2 < ν < 0.5. The upper limit case ν = 0.5 holds for isotropic, incompressible materials. Now, we can combine the material law (17) and the Green-Lagrange strain tensor E defined in Eq. (9) with the balance principles from above, leading to the following problem: Find displacement u such that holds, where Here, the boundary of our domain Ω is denoted by ∂Ω where homogeneous Neumann boundary conditions were assumed since no surface forces are acting on the boundary. Note that homogeneous Dirichlet values cannot be assumed anywhere on the boundary since it is a priori not known where the deformation will take place.
The latter is a considerable handicap, since the solution u is no longer unique, if we intend to omit the time derivative in the previous equation. Consequently, the time derivative from above has to remain as a stabilization which will require us to perform considerably more timesteps in practice.

Implementation of growth
As before, we use the deformation gradient decomposition approach (c.f., Eq. (1) main manuscript) to implement growth. Furthermore, we assume that the active deformation tensor F a may depend on local morphogen concentrations C, thus F a = F a (C), see also Eq. (6) main manuscript.
The key assumption behind the deformation gradient decomposition is that the intermediate state Ω a is stress free. Thus, any stress in the material body is generated by the elastic deformation F e . As a consequence, the elastic energy denisty Ψ and thus the first and second Piola-Kirchhoff stress tensors depend only on F e : and thus in contrast to the usual relation (16) and the material law (17). As a consequence, P e n a depicts the force in the current configuration Ω(t) given per unit surface area related to Ω a . However, the intermediate configuration is not a real domain, it might not even be continuous (e.g. if F a is only piecewise continuous) and is never physically attained. Furthermore, a surface area element in Ω a might not be uniquely defined.
Similarly, Σ e n a depicts the force in the intermediate configuration Ω a given per unit surface area in Ω a , which directly yields analogously to (16): Since Σ e is defined on the intermediate configuration and conservation of momentum cannot be written here, we intend to transform the balance of linear momentum equation to the reference domain Ω.
Transformation to the reference configuration yields the following relation between the second Piola-Kirchhoff stress tensor Σ and its elastic equivalent: Note that both arguments of the tensor Σ e have to be transformed from the intermediate configuration Ω a to the reference one. Thus, the Jacobian F −1 a of the transformation occurs twice. For a detailed derivation of this relation using tensor algebra we refer to the literature, for instance, see section 10 in Lubarda and Hoger [3].
We point out that the elastic Green-Lagrangian strain tensor E e and subsequently the Cauchy stress tensor σ e and the second Piola-Kirchhoff tensor remain symmetric, since Using the relation (22) we just derived, we can finally insert Σ e into our structural equation (19) which is given in the reference configuration. To derive our full mechanochemical system of model equations, we will now focus on the description of morphogen dynamics.

Full mechanochemical model equations
In Eulerian coordinates, the prototypic reaction diffusion equation for a concentration of signalling molecules c(x, t) reads with the velocity v(x, t): with the diffusion coefficient matrix D ∈ IR 3×3 , which we will further specify. The reaction term R(x, t) contains the mechanical feedback on the dynamics of signalling molecules and was presented in Eq. (4) (main manuscript). We transform Eq. (23) to the reference domain Ω using the deformation χ : Ω → Ω(t) and the transformations ∇ X C( of the gradient and the partial time derivative of a scalar-valued function C, respectively. The proof of these transformations is based on the chain rule and is given in standard textbooks, e.g. Richter [4]. Further, we employ the notation C = c(χ(X, t), t), v = v(χ(X, t), t) for the sake of a simple presentation.
The reaction-diffusion equation is now given in the reference configuration i.e. in the same framework as the structural equation (19) and reads As motivated in the main manuscript, we prescribe individual diffusion rates D N in normal (or radial), Lagrangian direction N = |X| −1 X of the tissue sphere and D T in tangential, Lagrangian directions T 1 and T 2 (c.f. Eq. (3) of the main manuscript). This yields the diffusion coefficient matrix is the transformation of the Cartesian coordinates to these systems and can be interpreted as a rotation. The matrices are given by However, before we combine both equations, we want to point out that we have so far neglected an important consequence of the deformation gradient decomposition as introduced above. Namely, it implies a separation of the time scale of growth / active deformations from the elastic one. Indeed, in biological tissues, growth takes place over a timescale of several hours or days whereas the elastic material response occurs within seconds. Obviously, we are interested in long-term active deformations (growth) rather than unrealistic, high frequent oscillations on the elastic time scale, which are also numerically costly to resolve i.e. it requires signigicantly smaller and thus considerably more timesteps. Consequently, we would like to neglect the time derivative ρ 0 ∂ tt u in the structural equation ("quasi-stationary formulation"). However, this leads to convergence problems of Newton's method, since the resulting problem is no longer unique as described above. Thus, we introduce the time derivative as a stabilization term with a stabilization parameter ∼ 0.1 to balance both effects.
Hence, the final coupled mechanochemical problem reads: Let Ω ⊂ IR 3 , be a (boundend) domain in two or three dimensions. Further, let the boundary of our domain be denoted by ∂Ω where homogeneous Neumann boundary conditions will be assumed. Then, find displacement u and morphogen concentration C, with initial conditions u(X, 0) = 0, C(X, 0) = C 0 such that in Ω, Again, µ and λ are the Lamé constants and D ∈ IR 3×3 is the diffusion coefficient matrix which allows for individual diffusion rates in radial and tangential directions as derived above. The reaction term R = R(Σ, E, F, C) contains the mechanical feedback on the dynamics of signalling molecules and was presented in Eq. (4) (main manuscript). For apical and basal constriction, the active deformation gradient is given by F a = Q TF a Q, wherê F a (C) depends on the local morphogen concentration and Q is a rotation matrix, see Eq. (6) of the main manuscript. As presented there, each biological cell constricts individually but in the same manner. The idea is now to prescribe the same active deformation for each cell and transform that tensor such that its orientation coincides with the orientation for the specific biological cell under consideration. Therefore, we have introduced local coordinate systemsX in the origin m of every biological cell such thatF a is identical in all these systems. The whole active deformation tensor F a is now given as follows: The a rotation Q i = Q| Ki to the systemX i of the considered biological cell K i , the application of the constriction F a (independent of i i.e. of the current cell) and the backwards rotation Q T i and the translation −Q T i m i of the argument. We have illustrated these mappings in Fig. A with the restriction to two dimensions to simplify the presentation. These mappings are easily extended to three dimensions defining Q as two rotations around coordinate axes.
Note that both rotations Q and R, introduced in this section, are very similar. On the one hand however, Q depends on the midpoint m i of the biological cell K i under consideration, i.e. Q| Ki = Q i (m i ), and is constant within cells. On the other, R(X) depends on each point X.

Derivation of mechanical invariants
A tensor A in three dimensions posseses three so-called principal scalar invarants I i , i = 1, 2, 3. The first invariant of a tensor A is always the trace and last invariant is its determinant. In three dimensions, this gives the tensor invariants where ijk denotes the Levi-Civita symbol. Consequently, if our mechanical feedback influencing the production of morphogen is merely depending on these invariants, it does not change with the rotation of the coordinate system. To be more precise, the coupling term will be objective for an objective tensor A. Obviously, it is crucial that physical quantities such as deformation, stress or strain are invariant relative to a change of observer which should not effect the properties of the material under consideration. For our purpose, it is thus reasonable to take the tensor invariants of these quantities as a mechanical feedback in our coupling R. For a better understanding of objectivity, we refer to the presention in section 5.2 in Holzapfel [1].
Clearly, the second Piola-Kirchhoff stress tensor Σ is intrinsically objective, since it is parameterized by material coordinates only. Furthermore, it can be shown that the first Piola-Kirchhoff stress tensor P and the deformation gradient F are objective, although they behave like vectors. However, one of their indices is described in the reference configuration, which makes them objective nonetheless. Simple argumentation shows that the Cauchy stress tensor σ and the Green-Lagrange strain tensor E are objective as well.
This gives us several potential candidates for the mechanical feedback. Since we are in the reference configuration, it is natural to choose the second Piola-Kirchhoff stress tensor Σ form the available stress measures. With regard to extensive numerical testing, we restrict ourselves to the determinant of the deformation gradient I 3 (F) = det(F) for the purpose of this publication, see Eq. (4) of the main manuscript. Note, that different feedback loops also showed promising results which are intended to be published in futur works.

Robustness of pattern formation
We have demonsrated that our simple, positive feedback loops lead to mechanochemical pattern formation (see Fig. 2 from the main manuscript for basal constriction and Fig. 3 main manuscript for apical constriction). In the following, we focus on the robustness of such patterns with respect to changes in the initial conditions, the model geometry and the model parameters/assumptions.

Robustness with regard to initial conditions
Firstly, we demonstrate, that our positive feedback loop produces the same number and quality of patterns if we change the initial morphogen distribution. Regarding pattern formation for basal constriction Here, morphogen is initially distributed with a spot at one side of the sphere. In the lower row, the 3D tissue body has been sliced just for the purpose of a better visualisation.
(c.f. Fig. 2 main manuscript) we now start with just one weak morphogen spot at the side of the sphere (instead of a stochastic morphogen distribution) and again obtain a stationary solution with the same number and size of mechanochemical patterns (Fig. B) where the solution still slightly changes after four days but is stationary after 25 days. Secondly, the same robustness with regard to a change in the initial conditions is observed for apical constriction: also here, starting with a weak morphogen spot at the side of the sphere (instead of a stochastic distribution) our feedback mechanism again results in the same number and quality of transient patterns and finally also in a gastrulation event with an annulus of high concentration around the invagination (please compare Fig. 3 from the main document and Fig. C).

Comparable patterns for apical and basal constricton
At fist sight, from Fig. 2-3 in the main manuscript, it might not be obvious that the mechanochemical patterns gained by our positive feedback loops for apical and basal constriction are actually strongly related: In Fig. D we compare simulation snapshots for both cases with an identical orientation of the coordinate system and at the identical time. Apparently, the curvature patterns are very similar where the chemical concentrations just are negatively correlated.
Yet, only for apical constriction it was possible to simulate a complete gastrulation event, c.f. Fig.  3 (main manuscript) and Fig. C. In general, these events seem to occur in both considered feedback loops since there is usually one dominating curvature/morphogen patch which evolve faster than others, for instance see Fig. 3 (main manuscript) after four days (apical constriction) or Fig. B after two days Here, morphogen is initially distributed with a spot at one side of the sphere. In the lower row, the 3D tissue body has been sliced just for the purpose of a better visualisation.
(basal constriction). Note that this statement neglects artificial configurations, e.g. of two identical, opposing initial morphogen circles for apical constriction and presumes that mechanochemical coupling is strong enough to produce large patterns. However, calculations for basal constriction break down at an earlier stage, i.e., Newton's method no longer converges, for instance, if we increase the morphogen production rate to obtain larger deformations and an eventual gastrulation rather than the stationary patterns observed in Fig. 2 (main manuscript). However, some of our calculations indicate that a mechanical feedback relying on stress might produce more stable gastrulation patterns (results not shown).
This gastrulation event for our feedback loop based on apical constriction means a significant topological change of the tissue geometry: eventually, we do not any longer observe a further increase in morphogen concentration but the invagination continues until the model breaks down (as discussed in the main manuscript). However, first simulations regarding an internal volume constraint are very promising: if we prescribe a surface force in normal, outward direction of the internal boundary of the tissue sphere which becomes stronger if internal volume is lost, this pressure stops the gastrulation event and seems to ensure a stationary result (c.f. Fig. F (c) ). Robustness regarding the model geometry and parameters As described above, changes in the initial condition lead to qualitatively and quantitatively comparable patterns. In contrast, changes in model parameters or geometry may also influence the number, regularity or size of patterns. One example is given by the tangential diffusion, which appears to smooth/regularise morphogen patterns: simulation results show that quartering the tangential diffusion leads to more but smaller patterns ( Fig. E (a)). In fact, tangential diffusion is not essential to produce patterns at all: Without any tangential diffusion we still obtain mechanochemical patterns, shown for basal constriction in Fig. F (a) (while this configuration obviously contains the modelling error that morphogen also no longer diffuses in tangential directions within biological cells). At first, these patterns are very tiny but as they deform surrounding tissue which is stretched, this stretch leads to morphogen production via our positive feedback loop and the patterns slowly extend (tangential diffusion also accelerates pattern formation).
We also investigated the influence of the model geometry on resulting patterns. Firstly, we considered a smaller tissue sphere consisting of only 384 biological cells (see Fig. E (b)). It appears that we obtain comparable patterns to those of our original domain in Fig. 2 (main manuscript) and Fig B, however there are less morphogen patches due to the smaller surface of this sphere.
Additionally, we observe that morphogen/curvature patches dominate larger parts of the tissue sphere if the tissue thickness is decreased (Fig E (c)). Contrary, there are fewer patches deforming larger parts of the tissue if the thickness of the sphere is increased, c.f. Fig E (d) (which in principle also holds for Fig E(b) where the relative tissue thickness is increased in comparison to our usual domain). In the "parameter setup" section of the main manuscript we stated that we essentially balance tangential diffusion and the morphogen production rate k 2 to obtain mechanochemical patterns. Whereas the effect of tangential diffusion was just discussed, the influence of the latter can be illustrated as follows. A larger production rate results in higher morphogen concentrations which in turn induce larger active cell-shape changes i.e. larger deformations. Analog to patterns in a thicker tissue sphere, these patterns dominate larger parts of the domain. For instance, we usually expect more patterns on the thinner domain in Fig E (c) which is outbalanced by the higher morphogen concentrations and thus larger deformations in this example.
In summary, we see, that smaller lateral diffusion and thinner domains or lower morphogen production result in more patterns; larger diffusion, thicker domains and higher morphogen levels lead in contrast to fewer and larger patterns.

Robustness regarding active deformations
Finally, we show that our simple feedback loops also lead to mechanochemical patterns if we change the nature of the active deformations. If we e.g. redefine the active deformation tensor such that it constricts one side of the cell without actively expanding the other side (i.e., the active deformation is not any longer volume preserving) we still obtain regular patterns (Fig. F (b)). Interestingly, morphogen patches are less spherical but rather develop to "double-patches", i.e., two superimposed patches constituting Here, we investigate the robustness of resulting patterns with respect to changes in lateral morphogen diffusion respectively non-volume preserving active deformations. (d) Simulation snapshots showing the gastrulation based apical constriction, which was stopped by an internal volume constraint. In the lower row, the 3D tissue body has been sliced just for the purpose of a better visualisation. together a dumbbell-shape.