A Mechanochemical Model for Embryonic Pattern Formation: Coupling Tissue Mechanics and Morphogen Expression

Motivated by recent experimental findings, we propose a novel mechanism of embryonic pattern formation based on coupling of tissue curvature with diffusive signaling by a chemical factor. We derive a new mathematical model using energy minimization approach and show that the model generates a variety of morphogen and curvature patterns agreeing with experimentally observed structures. The mechanism proposed transcends the classical Turing concept which requires interactions between two morphogens with a significantly different diffusivity. Our studies show how biomechanical forces may replace the elusive long-range inhibitor and lead to formation of stable spatially heterogeneous structures without existence of chemical prepatterns. We propose new experimental approaches to decisively test our central hypothesis that tissue curvature and morphogen expression are coupled in a positive feedback loop.


Introduction
During the development of biological tissues, different signaling molecules are responsible and indispensable for pattern formation and shape generation [1]. Since the seminal paper of Alan Turing [2] a variety of patterns in biological tissues have been studied using a framework of reaction-diffusion equations. These approaches assume that there exist diffusing signaling molecules, called morphogens, whose nonlinear interactions combined with different rates of diffusion lead to destabilization of a spatially homogeneous equilibrium and formation of spatially heterogeneous structures.
The idea looks counter-intuitive, since diffusion is expected to lead to a uniform distribution of molecules. Mathematical analysis of reaction-diffusion equations provides explanation of the phenomenon postulated by Turing. Patterns arise through a bifurcation, called diffusion-driven instability, in which a spatially homogeneous stationary solution looses stability for a certain range of diffusion coefficients and a stable spatially heterogenous stationary solution appear. The resulting structures can be monotone corresponding to the gradients in positional information or spatially periodic and their ultimate shape depends on the diffusive scaling, related to the size of domain.
The most famous embodiment of Turing's idea in a mathematical model of biological pattern formation is the activatorinhibitor model proposed by Gierer and Meinhardt [3]. One of the key ingredients of that model, responsible for Turing type dynamics, is that the inhibitor diffuses faster than the activator, i.e. the system is regulated by a short range activation and a long range inhibition [2,3]. However, in many developmental processes, dynamic and complex tissue topologies are likely to prevent the establishment of long range inhibitor gradients [4]. Furthermore, diffusion coefficients of typical morphogens are often found to be quite small [5], i.e. do not allow existence of significantly varying diffusion rates as required by the classical Turing mechanism. These observations suggest a search for a different inhibitory mechanism such as mechanical inhibition [4]. Moreover, mechanically based laws in morphogenesis appear to be promising and powerful alternatives to purely chemical models: The latter reduce macroscopic structures to 'blind' by-products of spatial chemical patterns and contradict certain experimental data [6,7].
The influence of morphogens on tissue mechanics (such as curvature) is well known. However, different studies show that the interplay is reciprocal and mechanical stress and physical forces (e.g. induced by tissue deformations) can also locally influence morphogen patterns and cell behavior [4,6,[8][9][10][11][12]. Furthermore, it appears that tissues may act differently depending on directions of applied stress [13].
Based on these observations, we propose a model for pattern formation in biological tissues, directly coupling the expression of a morphogen with tissue mechanics. Numerical simulations of the model based on our hypothesis reveal that N simple interplay between tissue mechanics and morphogen production can lead spontaneously to curvature and morpho-gen patterns with shapes depending on the size of the tissue; and N resulting patterns are insensitive to stochastic perturbations of initial conditions.
The proposed mechanism is a promising candidate to replace the ''missing'' long range inhibitor in the activator-inhibitor models. Due to its robustness, the ultimate shapes are reliably generated under various conditions. Therefore, the model can explain self-organization and de novo pattern formation in the systems with the initial state close to an equilibrium.
The mechanism is currently under experimental investigation. Furthermore, the presented mathematical model and numerical approach could be used in future work to investigate the interplay of different morphomechanical models such as those discussed in ref. [6].
Our mathematical model combines a reaction-diffusion equation for the morphogen with an elastic gradient flow for tissue mechanics. Expanding the ideas of Cummings [14] we assume that the morphogen locally induces positive curvature, and in turn, this curvature induces the expression of this morphogen (c.f. Fig. 1), i.e. that there exists a positive feedback loop. Above, the terms ''positive'' and ''negative'' curvature refer to outward and inward bending, respectively, compared to the initial curvature of a tissue. More accurately, experimental data indicate that it is not tissue curvature itself but curvature modulated tissue compression that influences gene expression [4,11]. However, in order to reduce the complexity of the model and to present the basic idea, we adhere in this paper to the simplified model. Models including tissue compressibility will be subject of future research.
Finite element simulations reveal that the postulated mechanism produces spontaneously a variety of curvature and morphogen patterns in an asymptotically stable way, i.e. insensitive to small changes in initial conditions. Based on a range of simulation results, we present detailed parameter studies of the model analyzing the rescaled parameter space. We identify parameters that control pattern-related scales, such as size, amount and curvature of appearing patches. Our results show that biomechanical interactions may constitute the missing link to the Turing long range inhibitor: positive curved domains induce negative curvatures at their edges; in these regions only morphogen degradation takes place (c.f. Fig. 1).
Furthermore, numerically obtained results show that dynamics as well as final curvature and morphogen patterns appear to be very similar to those observed during symmetry break in Hydra tissues. Thus, numerical results suggest that the presented mechanism could constitute a key mechanism for pattern formation in Hydra polyps. Further experimental work is needed to support the idea that the proposed mechanism indeed constitutes a key mechanism of pattern formation during embryogenesis. Corresponding experiments are suggested at the end of this paper indicated by simulations.

Model
LetX X : U?C be a time-dependent parametric representation of a closed cell sheet C, where the cell sheet is parameterized over the unit sphere adopting the ansatz of [15,16], i.e. U~S 1 . The morphogen level is described by the function w : U?R §0 ; for eachX X (u 1 ,u 2 )[C the concentration w(X X ) is identified with w(u 1 ,u 2 ). Hence, w is naturally moving with the deforming tissue, i.e. d t ½w~+ C ½w : L t ½X X zL t ½w, where + C is the surface gradient. For convenience of the reader, detailed definitions of the used geometrical quantities are given in the Table S1 available online.
To model the curvature-dependent elastic properties of a thin cell tissue, we use a modified Helfrich energy, i.e.
Here, H is the mean curvature and ds is the surface measure. k constitutes the bending rigidity reflecting the stiffness of the surface, and the spontaneous curvature H 0 reflects the preferred local surface curvature. The latter e.g. depends on the shape of biological cells within the tissue: If cells are symmetrical, the tissue energetically favors a flat geometry and H 0 vanishes. However, in the case of wedge-shaped cells, the tissue prefers to be locally curved, hence H 0 =0 holds. A detailed treatment considering the connection between the geometry of biological cells arranged in a tissue to Helfrich measures such as H can be found in [17]. Since laterally heterogeneous distributed molecules -such as morphogens-may influence the mechanical properties of the tissue, we assume (analogously to [18,19]) that the mechanical moduli k and H 0 are functions of the concentration w. In first approximation we take k(w)~azaw and H 0 (w)~bzbw as linear functions of w, where a §0. The original Helfrich energy F H~ð C k(H{H 0 ) 2 ds [20] with constant moduli k and H 0 has been widely used describing the mechanics of membranes. In the following, we assume a local tissue incompressibility. Adopting an energy point of view, the evolution of the tissue deformationX X up to time Tw0 in U|½0,T) is given by the following L 2 {gradient flow under the constraint of local incompressibility of the tissue: where d t is the total time derivative, L X is a kinetic coefficient, d dX X ½: denotes the variation with respect to the arbitrary vectorX X , and g is a local Lagrange multiplier [21] due to the constraint of incompressibility (g can be interpreted as pressure). Furthermore, ffiffiffi g p is the surface measure where g is the determinant of the first metric tensor (c.f. Table S1). Volume constraints are not considered (to account for the fact that in pattern formation experiments with Hydra tissues it appears that tissue spheres exchange internal fluid with their surroundings [22]). The gradient flow (2) leads to minimization of the free energy F tissue under the constraint of incompressibility (3). For detailed calculation of the variations of F tissue we refer to [19].
Let us now consider the dynamics of the morphogen w within the tissue. In contrast to the modeling of membrane dynamics showing fluid behavior with respect to lateral flows, the evolution of morphogens is modeled separately rather than obtained by a corresponding variation of free energy. Beside the basic assumptions concerning diffusion and degradation [5,23,24], we define the morphogen production as a function depending on the surface curvature. Furthermore, we assume that initially the tissue is arranged in a mechanically relaxed configuration with curvature H i~H (t~0,X X ). Experimental findings reveal that tissue deformations may induce morphogen production [4,8,9]. Thus, we postulate that the expression of w can be induced by local curvatures HwH i . If both, negative and positive curvatures had induced morphogen expression, i.e. H=H i , we would have expected oscillations instead of stable patterns. Hence, w induces local positive curvatures via H 0 (w) and positive curvatures induce locally the expression of w via appropriate reaction-terms, constituting a positive feedback loop (c.f. Fig. 1). Such positive Figure 2. A-B: Simulations of spontaneous tissue pattern formation starting from stochastically distributed morphogen on a sphere at different time steps. A: For a strong coupling between curvature and morphogen expression the system shows strong budding (numerics do not perform properly for tw20 since the approach is only able to handle medium geometry changes; the likely equilibrium shape consists of fully budded morphogen patches). B: For a weak coupling between curvature and morphogen expression the system reaches a symmetric mechanical equilibrium at tw50. C-D: Equilibrium patterns are insensitive to different choices of initial conditions. C: Plot of total morphogen level ð feedback loops between chemical and mechanical processes during development have just been reported [12].
Using Michaelis-Menten kinetics (assuming the existence of a maximal expression rate of the w-promoter) and defining H §0 :~maxf(H{H i ),0g, we obtain the dynamical equation for w: with constants C,d,f,v [ R §0 . Hence, the model is given by a nonlinear PDE system of fourth order, coupling the gradient flow for tissue mechanics (2) under the constraint of local area incompressibility (3) with the reaction-diffusion equation (4) for morphogen dynamics.

Nondimensionalization
In the following subsection the free energy (1) and corresponding dynamic equations (2)-(4) are nondimensionalized. This reformulation allows us, on the one hand, to identify essential parameters and characteristic properties of our system. On the other hand, it enables us to study the relative dependence of different parameters on spatial and temporal scales. To nondimensionalize the model, we replace each variable v with a dimensionless quantity v c which is scaled with a characteristic unit of measure E v .
Free Energy: We setX X~E xX X c with ½X X ~½E x ~m which implies as well as k(w)~E k k c (w) with ½E k ~½k(w)~k B T. It follows: A consistency check shows that ½F has units of energy as expected, i.e. F~E k F c with ½F c ~1. Calculating the variation of F we obtain for the strong formulation with ½ dF

Surface dynamics:
We set t~E t t c with ½t~s. Thus it holds Morphogen dynamics: Choosing E t and E x appropriately we can always guaranteẽ d d~R~1. Furthermore, simulations show that equilibirum patterns are independent of E w (c.f. Fig. 1 C-D). Hence, choosing E w appropriately we can always yieldb b~1. Moreover, restricting ourselves to the case k:1 (i.e. a~1 and a~0) andṽ v~1 as well as choosing b~H i for biological reasons the total parameter space can be reduced to three independent constants, namelyC C,f f and L Xc . Finite element approximation and parameter setup Using the outlined modeling approach, we investigate numerically the effect of different model parameters, system sizes and initial conditions on pattern formation. To do so, we approximate equations (2)-(4), and more accurately, the nondimensionalized form derived in the previous subsection, to study systematically spontaneous pattern formation induced by the proposed mechanism. The corresponding total parameter space can be reduced to three independent constantsC C,f f and L Xc .
For spatial discretization we use a biqubic mixed parametric finite element method. The tissue surface is discretized using a quadrangular grid with 1664 grid points. For time discretization we use an adaptive semi-implicit Euler scheme. For further details to the computation scheme we refer to ref. 16, 19. Related to the early shape of an embryo (a blastula) we assume as the initial tissue geometry the sphere S R (0) with R~1, choosing an appropriate nondimensionalization, c.f. previous subsection. Moreover, we use a stochastic distribution for the initial morphogen concentration. The corresponding values at each discretization point are chosen using the standard random generator provided by C++. w(t~0) is uniformly distributed in an interval ½0,z with the average Sw 0 T~ð In the presented study, we restrict our studies to the case of morphogens influencing tissue curvature but not tissue rigidity, setting k:1 (i. e. a~1 and a~0). Moreover, we setṽ v~1, since neglecting long time behavior, i.e. nonlinear effects, it is sufficient to vary b=v and b~H i . This latter represents the assumption that the tissue is initially arranged in a mechanically relaxed configuration. Furthermore, due to an appropriate nondimensionalization, we can always guaranteed d~b b~1 (c.f. previous subsection).

Results and Discussion
Here, we will work only with the nondimensionalized form of equations (2)-(4) (c.f. previous sections). To investigate the variety of patterns produced by the presented model, we have performed more than 250 simulations using the nondimensionalized param-etersC C,f f and L Xc (setting a~ṽ v~b b~R~1, b~H i and a~0, c.f. previous sections). Furthermore, we have considered different values of Sw 0 T. Depending on the exact choice of parameters, simulation results show two different types of long time behavior: either the system tends to equilibrium with a symmetric medium curved pattern and seems to be stable at least for simulated times (c.f. Fig. 2 B), or strong budding appears (c.f. Fig. 2 A) and persists as long as the numerical scheme functions (the adopted numerical discretization is only able to handle medium geometry changes); it is likely that the limiting shape consists of fully budded patches, as experimentally observed e.g. during the formation of new body axes in Hydra polyps. For the corresponding simulations we have set L Xc~C C~1 as well asf f~4 in C butf f~8 in B. Using different initial conditions (i.e. different values of Sw 0 T=0 on an undeformed sphere) we find that the equilibrium pattern is independent of the exact choice of initial conditions (c.f. Fig. 2 C) and corresponding equilibrium patterns are very similar (c.f. Fig. 2  D). Even starting with no morphogen (Sw 0 T~0) and stochastically perturbed geometry yields the same result (not shown). Hence, the presented mechanism appears to be asymptotically stable and does not require the existence of any prepattern.
Checking systematically the influence of the varying model parametersC C,f f and L Xc on the corresponding emerging patterns, reveals the following relationships: N The distance between the patches (and hence the number as well, but not the size) and their curvature can be controlled by the strength of coupling between curvature and morphogen expression, i.e. by the curvature dependent productionf f: In Fig. 3 B the results of a stepwise increase fromf f~2 up tof f~9 keeping R~1,C C~0:01 and L Xc~0 :0001 are shown.
N Changing L Xc does not result in striking differences in pattern formation but only influences relevant time scales, i.e. the choice of the observed steady state is independent of L Xc . This is expected since possible stationary states are independent of L Xc (the whole right hand side of equation (2) is multiplied by L X ). Having in mind that the approach here is based on the assumption of an overdamped motion, L X can be interpreted in terms of viscosities. That is, patterns are independent of the specific viscous damping.
Investigating the dynamical behavior of the simulations reveals that all performed simulations typically follow three qualitatively different paths during evolution (c.f. Fig. 2 A-B; for the corresponding simulation we have setC C~0:01,f f~8, and L Xc~0 :01.): First, smoothing of initially stochastically distributed morphogen distribution; second, appearance of new and larger slightly curved morphogen patches; and third, distinct curvature patterns coinciding with high morphogen levels, either stabilizing at a symmetric pattern or strongly budding. Interestingly, microscopy data of evolving Hydra-reaggregates show similarities to these dynamical patterns [25].
Apart form the discussion concerning existence and nature of the morphogens (see Introduction), there exist previous mechanochemical models [26,27] conceived to explain budding processes in polyps. Also, as it is known, reaction-diffusion systems exhibit similar qualitative behavior, although under differing assumptions. In order to motivate experiments proving (or disproving) the proposed mechanism, we have performed 'virtual experiments' shown in Fig. 4. Since it is already known that morphogens are able to influence tissue curvature, we focus in these experiments on the effect that tissue curvature induces morphogen production (depending on the direction of tissue curvature): Forcing the tissue to bend locally outwards (e.g. using micropipettes) activates the morphogen level strikingly (Fig. 4 B) whereas inward bending prevents the establishment of morphogen patches in these regions (Fig. 4 C). However, no applied force yields randomly distributed morphogen patches (Fig. 4 A). In all corresponding simulations we have setC C~0:01,f f~8, and L Xc~0 :01 as well as circular pulling and pressing force-terms in B and C, respectively.
To summarize, in this contribution, we have proposed a new non-Turing type model for early pattern formation in tissue development. Based on recent experimental findings, the key assumption is a positive feedback loop between tissue curvature and morphogen production. We have shown numerically that this simple mechanism itself leads to various morphogen and curvature patterns resembling those observed in experiments involving spherical aggregates of Hydra-tissue. One of the aims of this paper is to motivate further experimental research to validate the presented mechanism for early pattern formation.
If it is validated by direct experimentation, for example along the lines suggested in Fig. 4, the proposed mechanism will constitute an essential step in the evolution of an initial homogeneous tissue sphere to a complex organism -one of the greatest current mysteries in biology [28,29].

Supporting Information
Table S1 General mathematical notations and definitions used in this paper. (PDF)