Beyond Turing: mechanochemical pattern formation in biological tissues
 Moritz Mercker^{1}Email author,
 Felix Brinkmann^{1, 2},
 Anna MarciniakCzochra^{1} and
 Thomas Richter^{2}
DOI: 10.1186/s1306201601247
© Mercker et al. 2016
Received: 18 December 2015
Accepted: 20 April 2016
Published: 4 May 2016
Abstract
Background
During embryogenesis, chemical (morphogen) and mechanical patterns develop within tissues in a selforganized way. More than 60 years ago, Turing proposed his famous reactiondiffusion model for such processes, assuming chemical interactions as the main driving force in tissue patterning. However, experimental identification of corresponding molecular candidates is still incomplete. Recent results suggest that beside morphogens, also tissue mechanics play a significant role in these patterning processes.
Results
Combining continuous finite strain with discrete cellular tissue models, we present and numerically investigate mechanochemical processes, in which morphogen dynamics and tissue mechanics are coupled by feedback loops. We consider three different mechanical cues involved in such feedbacks: strain, stress, and compression. Based on experimental results, for each case, we present a feedback loop spontaneously creating robust mechanochemical patterns. In contrast to Turingtype models, simple mechanochemical interaction terms are sufficient to create de novo patterns.
Conclusions
Our results emphasize mechanochemical processes as possible candidates controlling different steps of embryogenesis. To motivate further experimental research discovering related mechanisms in living tissues, we also present predictive in silicio experiments.
Reviewers
Reviewer 1  Marek Kimmel; Reviewer 2  Konstantin Doubrovinski (nominated by Ned Wingreen); Reviewer 3  Jun Allard (nominated by William Hlavacek).
Keywords
Morphogens Tissue morphogenesis Development Pattern formation Mechanochemistry Tissue mechanics Mechanotransduction Reactiondiffusion Longrange inhibitionBackground
with zeroflux boundary conditions. Variables c and h denote concentrations of two morphogens, called activator and inhibitor, respectively. The parameters ρ _{ c } and ρ _{ h } are production/removal rates, and d _{ c } and d _{ h } denote the diffusion coefficients.
Mathematical analysis of the equations provides explanation for the phenomenon postulated by Turing, for details see [33, 39]. The mechanism is related to a local behavior of solutions of a reactiondiffusion system in the neighborhood of a constant stationary solution that is destabilized through diffusion. Patterns arise through a bifurcation, called diffusiondriven instability (DDI). The specific choice of model kinetics and distinct different diffusion rates d _{ h }>>d _{ c } are necessary to induce the DDI (c.f., Fig. 5 a). Consequently, after more than 60 years of research, molecular identification of the Turingtype molecules during morphological events remains an exception rather than a rule [16, 26]. Especially, appropriate candidates for fast diffusing long range inhibitors are usually missing. Interestingly, recent results suggest that in general the principle of “longrange inhibition/short range activation” holds. However, specific realizations seem to be versatile: e.g., primarily mechanical mechanisms have been shown to produce periodic tissue patterns as well (reviewed in [16]). Mechanical cues such as traction, rigidity, or compression underlie functions in patterning [7, 24, 40, 51], which have been ascribed to diffusing morphogens within Turingtype models. The active role of tissue mechanics in tissue patterning becomes justified by the fact that an increasing amount of data shows how gene expression and signaling cascades sense and process mechanical cues (for reviews, c.f. [5, 6, 12, 19, 30–32, 41, 44, 56, 57]). Taken together, these data suggest that a close interplay between chemical and mechanical processes within tissues may lead to pattern formation in many cases.
Despite the large number of experimental data, modeling approaches explaining de novo pattern formation as results of mechanochemical processes are still rare. Existing models integrate the tissue mechanics only in a simplified way [36, 37, 52]. Thus, the need for new modeling approaches has been recently stressed [18, 55]. A general challenge in such mechanochemical models is the description of the mechanical part of possible mechanisms, since it is nonintuitive from a modeling/numerical point of view. Already Turing suggested in his seminal paper the consideration of mechanical aspects in pattern formation, but restricted his own study to the chemical processes, since ”...the interdependence of the chemical and mechanical data adds enormously to the difficulty” [53]. Currently, modeling and computation approaches describing mechanical aspects of morphogenesis have reached an advanced level (for reviews, c.f. [48, 58]), making the new models feasible.
In the current study, we propose and investigate integrated mechanochemical models to investigate tissue pattern formation. The models are based on the finite strain theory [17] and are extensions of the works as presented in [1]. Based on experimental data, we model different mechanochemical feedback loops by coupling diffusing morphogens with different mechanical cues, such as compression/stretch [5], strain/cellshape [29, 56] or stress [28, 41, 42]. We numerically investigate resulting equilibrium patterns and their dependence on diffusion rates, initial conditions or the system size. Finally, we propose experimental approaches to decisively test analogous mechanochemical interplays in vivo.
Results and discussion
In this section, we present and discuss several simulation results which are based on the proposed mechanochemical tissue models. In particularly, we combine continuumbased finite strain models with discrete models for biological cells. Furthermore, the simulated tissue geometry is based on an annulus geometry, so that we consider a thin loop of tissue confined to the plane. For more details regarding the model geometry, the modeling and simulation approach, the underlying biological/biophysical assumptions, as well as the detailed parameter setup, we refer to the “Methods” section.
Mechanochemical pattern formation
We want to point out that the ability to form patterns in all three mechanochemical feedback loops does not critically depend on the initial conditions: Each presented feedback loop yields spontaneous emergence of stable patterns starting from a stochastically as well as a gradually distributed morphogen (results not shown).
In all presented simulations, we consider an apical or a basal constriction as an active part of the deformation, since such deformation processes appear to be frequently involved in tissue morphogenesis [34, 38, 50]. However, other possible active processes – such as isotropic or anisotropic tissue growth – can also play an important role in morphogenesis, especially in embryogenesis [2, 10], and are not considered in the models presented. Although detailed results are postponed to future works, our first simulation studies indicate that also here, different mechanochemical feedback loops spontaneously create robust patterns (results not shown). Thus, mechanochemical pattern formation does not critically depend on the exact nature of the active deformation.
Our simulation results furthermore reveal that the resulting equilibrium patterns can qualitatively differ, depending among others on the specific type of the mechanical cue involved in the feedback. Most frequent patterns show at least two morphogen/deformation patches arranged in a symmetric manner, such as footsole shapes (Fig. 1 a, right) or ellipses (Fig. 1 b, right), which can be generated by each of the three feedback loops. But also patterns with only one morphogen/curvature patch are possible (Fig. 1 c, right), resembling those appearing during gastrulation event in embryogenesis [23]. Interestingly, the latter geometry we only observed in the case of the stressmediated feedback loop, and starting with a gradually distributed morphogen. Thus, the symmetry and type of final patterns may critically depend on initial conditions and the type of the mechanochemical feedback loop.
It is worth mentioning that the patterns resulting from simulation of tissues restricted to the 2D plane may distinctly differ from those of full 3D tissues, for several reasons: From the mechanical point of view, we assume that the 3D nature of a blastula coupling cells and introduction of additional curvature may affect model outcomes. From the chemical point of view, a process such as circumferential pattern formation lacks one dimension, which implies that no distinction is made between stripes and spots. Thus, the models investigated in this work can be treated as test cases showing that a variety of stable patterns is produced from mechanochemical feedback loops, while for a comparison with experimentally observed structures 3D models are needed. We postpone the extension to the 3D models to future research.
Mechanics versus diffusion in pattern formation
We emphasize the crucial role of tissue mechanics in shortrange activation and longrange inhibition in all feedback loops considered: If cells are locally constricted due to increased local morphogen abundance, two simultaneous events promote subsequent pattern formation, 1) the local deformations induce an increased morphogen production at the same place, which leads to even stronger deformations and as a result we obtain a shortrange activation of the morphogen; 2) the local deformations immediately lead to stress/strain/compression at surrounding tissue areas, which inhibits the local production of the morphogen at these places. As a result, in this latter case we obtain a mechanically mediated longrange inhibition, and local high morphogen concentration patches can inhibit other patches over long ranges.
Predictive experiments
 1.
the existence of a morphogen locally changing tissue mechanics, as well as
 2.
a reverse signal from mechanics back to chemistry, e.g., a mechanical cue changing effective morphogen spread and/or production.
Having selected a certain morphogen candidate, each causality direction can be subsequently separately assessed: To prove (1), the morphogen of interest could be locally and ectopically overexpressed, with a subsequent search for the induction of colocalized mechanical patterns (such as deformations, c.f. [45]). Although the visualization of mechanical patterns different from deformations (such as stress or compression) is still challenging (e.g., [42]), recent technical developments are promising [22, 48]. Analogously, to assess item (2), the tissue can be actively deformed [11, 48], combined with a subsequent screen for accordingly aligned morphogen expression patterns.
Conclusions
To summarize, in this publication we have presented and investigated a class of models of pattern formation in biological tissues, namely mechanochemical models. We simulated three different simple feedback loops between mechanical cues and morphogen dynamics. In all three cases, we have shown that these type of interactions lead to a de novo formation of stable chemical and mechanical patterns. These results indicate that there are various possible ways in which a tissue may produce patterns by simple mechanochemical interplays. Furthermore, several of our predictions from simulations appear to be robust regarding changes in the mathematical model, such as changes in the type of the mechanochemical feedback loop (i.e., changes in the involved mechanical invariant or active deformations process), or in the exact choice of initial conditions.

In contrast to the Turingtype models, within presented feedbackloops, simple (linear) interaction terms and moderate diffusion rates are sufficient to produce robust de novo patterns. This makes the evolution of such mechanisms more likely.

Within presented mechanisms, mainly mechanical cues undertake functions in shortrange activation and longrange inhibition, chemicals (morphogens) play only a mediating and thus subordinate role. This could explain why the experimental identification of relevant molecular candidates is still missing in many cases.

Especially in developmental steps where tissue deformations play an important role (e.g., during gastrulation), presented feedback loops display a natural and simple way for the tissue to control the progress and success of deformations, since the latter directly retroact on morphogen dynamics.
One of the main aims of this paper is to motivate further experimental research in order to validate mechanochemical mechanisms for tissue pattern formation. Recent experimental observations [4], theoretical results [37], as well as technical developments [22, 48] are however promising. If mechanochemical mechanisms are validated in tissues by direct experimentation (for example similar to the methods as suggested in Fig. 3) this will constitute an essential step in the understanding of embryogenesis  one of the greatest current mysteries in biology [23, 39].
Methods
Model geometry
Combining continuous and discrete tissue models
The presented modeling approach combines continuous (finite strain) modeling techniques with a discrete model for the shape and position of biological cells. The explicit distinction of biological cells from numerical cells (finite elements) is necessary to appropriately describe active shape changes of biological cells, without being restricted by the local resolution of continuous mechanical processes determined by finite elements.
As an example, the active cell deformation processes of apical or basal constriction (c.f., following subsections) cannot be described within a purely continuous framework: although applied to the overall body, active deformations (given by F _{ a }, c.f. following subsections) here also depend on the discrete model of biological cells in the sense that the direction of deformation jumps at these boundaries (c.f., Fig. 6 a). These discontinuities represent the fact that the cytoskeleton is pulling from both directions at the boundary regions between biological cells. However, the overall deformation, F=F _{ e } F _{ a }, includes also the simultaneous passive (elastic) response, maintaining (among others) the continuity of the tissue. In terms of biophysics, the discrete part of our model represents the cytoskeleton associated with the plasma membrane of each biological cell (actomyosin cortex), where active deformations may be driven by myosin motor proteins. In contrast, the continuous model part represents the body of the cell lumina.
Deformation gradient decomposition and model equations
where F _{ a }(c) is a prescribed active deformation depending on the morphogen concentration c and F _{ e } is the passive elastic response to these active changes, which ensures the continuity of the overall deformation F and minimizes the mechanical stress (see also Fig. 6 a). Numerically, this means to replace F _{ e } by \(\mathbf {F} \mathbf {F}_{a}^{1}\) and to solve for the overall deformation F, applied to the continuous tissue body. Thus, although F _{ a } may show discontinuities at boundary regions between biological cells, the smoothness of the overall deformation is naturally given by the interplay between active and passive processes leading to the smooth finite element solution determined by F.
Furthermore, the active deformation F _{ a } (influencing the shape of the biological cells, c.f. Fig. 6 a) is directly determined by morphogen concentrations c. This is a common assumption in the mechanochemial modeling of tissue morphogenesis [1, 36, 43], and is based on the idea that morphogens locally induce remodeling of the cytoskeleton, which again determines the cell shape. In accordance, various experimental results show that different active deformation processes (such as cellshape changes or local tissue growth) can be directly induced by expression of signaling molecules [9, 27, 35, 54]. However, the final tissue shape is not entirely determined by F _{ a }(c) but rather by the interplay between F _{ a }(c) and the elastic response F _{ e } (c.f., Fig. 6 a). In this way, local morphogen concentrations (respectively the resulting F _{ a }(c)) can induce complex patterns of stress, strain, and compression within the surrounding tissue (Fig. 6 b).
with the Lamé constants μ and λ. The reaction term R=R(Σ _{ e },E _{ e },F _{ e },c) will incorporate the feedback of the mechanics on the morphogen level. With the exception of the second predictive in silicio experiment which is presented further on, the right hand side is set to f=0, i.e. external forces are usually not considered.
We emphasize that we are interested in a quasistationary state where the time derivative in the structural equation is not considered. In practice however, the time derivative multiplied by ε=0.1 is used for stabilization of the numerical scheme. The specific choice of ε does not significantly change the final pattern. Furthermore, the evolution equation with nonzero time derivative yields uniqueness of the numerical solutions, since homogeneous Neumann boundary conditions are assumed (for more details, we refer to the section “Finite element (FEM) approximation”). Consequently, the solution of the quasistationary system is unique up to translations and rigid body rotations. Such solutions might affect the convergence of the numerical scheme. Therefore, we consider a small parameter ε to balance the oscillatory behavior introduced by the time derivative with rigid body rotations experienced in the quasistationary setup. With such stabilization term, oscillations are largely overdamped with a negligible loss of convergence speed.
Mechanochemical events
Active deformations Mechanochemical pattern formation requires an interplay between active and passive chemical and mechanical processes. Especially, the morphogen concentration c couples into the structural equation through the active part of the decomposed deformation gradient (Eq. (1)).
where k is a constant and \(\left (\hat {\mathbf {X}}_{0},\hat {\mathbf {X}}_{1}\right)\) are the 2D coordinates in the cellwise reference system. For positive values of k, this results in apical constriction and for negative values in a basal one.
where the last integral vanishes since integration over x _{1} cancels out.
Mechanotransduction The role of mechanosensitive mechanisms controlling chemical cellular processes has been extensively studied within the last decade [12, 30]. Different types of mechanical cues have been shown to influence gene expression (such as morphogen production), namely stress [28, 41, 42], compression/stretch [5] but also geometrical constraints determining the strain/cellshape [29, 56]. In terms of continuum mechanics, these three cues can be expressed via invariants of the corresponding tensors. Tensor invariants were chosen since their values do not change with the rotation of the coordinate system, which is equivalent to a rotation of the initial morphogen concentration. Thus, using tensor invariants, the feedback and consequently the solution of our system will rotate accordingly to how the initial conditions were rotated.
Since there are two tensor invariants available in 2D, there are a total of six possible candidates resulting from mechanical stress, compression/stretch and strain.
 1.
the determinant of the deformation gradient I= det(F), which has the physical interpretation of compression or stretch. More precisely, \(\det (\mathbf {F}) = \frac {V_{a}}{V_{0}}\) is the ratio of the deformed to the initial volume;
 2.
the isotropic strain, which is the trace of the GreenLagrangian strain tensor, i.e., I= tr(E). This measure represents the hydrostatic strain which is the displacement between particles in the principal coordinate direction inside the tissue relative to a reference length; and
 3.
the determinant of the elastic Cauchy tress tensor \(I = \det (\mathbf {\sigma }_{e}) = \det \left (J_{e}^{1}\mathbf {F}_{e} \mathbf {\Sigma }_{e} \mathbf {F}_{e}^{T}\right)\). Let us note that σ _{ e } is used since the determinant of the antisymmetric tensor Σ _{ e } is not an invariant. σ _{ e } expresses the stress which biological cells are experiencing as a result of their active deformation.
Conceptually, stress is an internal force acting on a boundary per unit area of this boundary, while strain or compression are measures of deformation. For an illustration of the qualitative differences between these three different tensor invariants, we refer to Fig. 6 b. It appears that differences between compression and the strain invariant are usually small and are apparent only at the boundaries of biological cells.
with positive constants k _{1},k _{2},k _{ m }>0. Here, k _{1} is equal for all three considered feedback loops and causes a constant morphogen degradation in the entire tissue. The first term on the right hand side of Eq. (7) represents the impact of mechanics on the morphogen production: If the value of the respective tensor invariant I is positive, morphogen is produced. Else, no feedback is generated. For large values, this term converges towards k _{2} which can consequently be interpreted as the maximal morphogen production rate. k _{ m } is the Michaelis constant and is the value of I where half of the maximal production rate k _{2} is reached.
MichaelisMenten kinetics has been used to reflect the fact that there exists a maximal expression rate of the morphogen promoter. Thus, morphogen production saturates for large values of I, i.e., for large deformations, stresses, or strains.

local morphogen accumulation leads to local tissue deformations (in terms of apical or basal constriction), and

different mechanical cues (such as compression, strain, or stress) can induce the production of this morphogen in turn.
Predictive in silicio experiments
To demonstrate a possible experimental outcome indicating the existence of a mechanochemical feedback within a tissue, we have performed two predictive in silicio experiments using the strainmediated feedback loop as an example. For the detailed motivation of the following two experiments we kindly refer the reader to the Results and discussion section.
with k _{ x }=9·10^{9}. The first component f _{0} of this volume force is always zero. Since the center of the tissue is in the origin, the second component, f _{1}, constitutes a pull in negative x _{1} direction, i.e., downwards. The pull is localized around x _{0}=0 where it is strongest and decreases exponentially for larger absolute values of x _{0}; see also Fig. 3 c. Note that x _{1}<0 excludes a second pull downwards on the top of the tissue. Also, it is strongest in at the beginning of the deformation process and decreases exponentially in time.
Finite element (FEM) approximation
For discretization in space, different feedback loops required different levels of meshrefinement. On the one hand, for the compression and strainmediated ones, the entire domain Ω has been split into 4096 finite elements, see Fig. 4 a. Here, each biological cell is thus represented by 64 finite elements. Especially, isoparametric Q1 finite elements have been used (where Q1 depicts an approximation with polynomial ansatz functions of degree one).
On the other hand, the stressmediated feedback loop depends not only on the smooth solution u but also on the nonsmooth measure \(\mathbf {F}_{a}^{1}\), which required a mesh of at least 262144 cells in order to properly resolve the stressmediated feedback within the Q1 ansatz. It turned out that shorter computation times were possible using Q2 finite elements instead, i.e., with polynomial ansatz functions of degree two. This approach, combined with local mesh refinement, reduced the total amount of cells to a total of 17944.
Parameter setup
For the following calculations we have used E=100P a and ν=0.4, similar to the assumptions made in [8], and ρ ^{0}=1000k g m ^{−2} for the initial mass distribution. Furthermore, we have always set: D∼10^{−14} m ^{2} s ^{−1} for the diffusion coefficient, k _{1}∼10^{−4} s ^{−1} for the degradation rate of the morphogen level in the entire domain and k _{2}∼10^{7} m o l m ^{−2} s ^{−1} for the maximal morphogen production rate. The Michaelis constant was set to k _{ m }=0.1 which means that for positive feedback I=k _{ m } half of the production rate k _{2} is reached.
Finally, uniformly distributed random concentrations for each biological cells or a morphogen gradient were used as initial conditions. In both cases, the morphogen concentration was prescribed in the interval c∈[0,10^{9}]m o l m ^{−2}. For visualization in the following results, initial conditions were transformed into the interval c∈[0,1]m o l m ^{−2}. The scale of the morphogen does not influence results, since only the constant k which determines how strong the morphogen concentration couples into the active deformation gradient (see Eq. (5)), has to scale in the same manner. It was set to k∼10^{−6} m o l ^{−1} m.
Author’s response
Revision 1
Reviewer 1: Marek Kimmel – Rice University, Houston, Texas
Summary:
The manuscript presents a novel explanation of the morphogen in pattern formation mechanisms in development. It is proposed that since the putative purely chemical morphogens have been difficult to find, it is necessary to look for other possibilities. One such possibility is offered by mechanochemistry. Authors carefully demonstrate, using a mathematical model, how tissue mechanics provides a signaling modality. This is a clearly written and interesting paper.
 1.
I like the way the paper is written. The only suggestion I may have is that since the readership of Biology Direct may be somewhat less familiar with specialized models of pattern formation, the authors provide a descriptive, mathematical, or graphical representation of the "classical" Gierer and Meinhardt (or similar) model, so that the failure to find the putative morphogen is explained using a more specific example.
We have added a more detailed mathematical description of the classical GiererMeinhardt model as well as corresponding equations to the manuscript (page 2, paragraph 12). Furthermore, we have added a a graphical representation of these models to Fig. 5 (c.f., Fig. 5 a ).
Reviewer 2: Konstantin Doubrovinski (nominated by Ned Wingreen) – Universität des Saarlandes, Saarbrücken, Germany
Summary:
In their manuscript the authors present a mathematical model of morphogenesis that involves mechanosensitivity, i.e. a model where the dynamics of morphogen gradients depends on stress distribution in the tissue. It has long been known that biochemical processes can influence tissue mechanics, for example through control of molecular motors. In the recent years it is becoming increasingly realized that mechanical stresses can also influence biochemistry. In particular, it has been experimentally demonstrated that cell division rate and planar cell polarity in certain tissues can be influenced by external forces applied to those tissues. The model proposed by the authors is concisely summarized on page 8. The corresponding equations express conservation of momentum in the tissue and conservation of mass of the morphogen. Morphogen concentration drives tissue deformation through a contribution to the active deformation gradient (Eq. 3). Mechanosensitivity enters the description through a source term in the equation that describes morphogen dynamics. This source is assumed to be a function of some invariant of stress (or strain) tensor. The authors consider three models that incorporate three possible scalar invariants that might influence morphogen dynamics (listed on page 9). The paper discusses a subject of significant interest to biology community and in my opinion is appropriate for the intended readership.
On the technical side, there are a number of points that I think require clarification.
 1.
In the equation proposed by the authors the active deformation is directly determined by morphogen concentration. To me it appears more natural to assume that stress (rather than deformation) is a function of concentration. Could the authors comment on this?
We have added a corresponding paragraph including a discussion of this topic as well as additional references to the manuscript (c.f., page 9).
 2.
It seems to me that active deformation (denoted by F _{ a }) is a tensor and should thus transform like one (just like scalar invariants describing the coupling of stress and morphogen concentration dynamics transform like scalars). Equation 3, however, does not appear an isotropic tensor to me. I would rather expect this tensor to be a scalar multiple of identity matrix, where the scalar factor may depend on morphogen concentration. In any case, stress must be a symmetric tensor. The authors must explicitly show that their particular choice of F _{ a } ensures that stress transforms as a tensor and is symmetric.
The deformation gradient F _{ a } itself must not be symmetric or isotropic, it is just the gradient of any active mapping. Considering isotropic growth, F _{ a }=c I would indeed be the correct choice. We have added on page 11 a section explaining the setup of F _{ a } in the case of basal or apical constriction in detail. Further, we have shown, that the elastic GreenLagrangian strain tensor E in Eq. (4) is always symmetric, c.f. page 11. Hence, every stress tensor derived from E will be symmetrix.
 3.
In an actual biological setting I would expect the dynamics to be strongly overdamped implying that the dv/dt term in momentum balance equation must vanish. Was this true of the simulated regime?
We consider a quasistationary system by considering dv/dt only as a stabilization term, providing uniqueness of solutions. We have clarified this by adding a corresponding paragraph to the manuscript (c.f., page 10, last paragraph, as well as Eq. (3)). Actual biological settings will be damped by various factors, that are not present in our simplified models.
 4.
The section on the details of insilicio simulation is much too short to fully appreciate the details of the implementation. I think this section must be substantially extended.
As suggested by the referee, we have substantially extended this section (c.f., page 1314).
Reviewer 3: Jun Allard (nominated by William Hlavacek) – University of California, Irvine
Summary:

clarity in describing their model, especially its geometry, and

the generality of their results. In addition,

the conceptual distinction between the three mechanical invariants should also be clarified.
 1.
The authors should edit the text to clarify that this is a thin loop of tissue confined to the plane. This is important because the restricted geometry limits the significance of the results.  Add a paragraph at the beginning of “Results and discussion” stating the model geometry. Two or three sentences would be sufficient.  Delete phrases like “tissue sphere” (p.3), which risk being misleading.  Change the dimensionality language, for example at the top of p.4. The tissue is quasi1D in a 2D domain, so describing it as 2D is ambiguous. The sentence “in 2D domains there is no difference between stripes and spots” should be changed to something like, “In a thin strip, there is no difference between stripes and spots”. Figure 4 also risks being misleading, since the model is strictly continuumbased, and at no point are “biological cells” employed in the model. Indeed, Fig. 4 could be removed completely.
As suggested by the referee, we have revised the entire manuscript regarding the dimensionality language (see for example page 5, second paragraph), we have added additional information to the beginning of the “ Results and discussion ” section (page 3, last paragraph), and have added information to the “ Model geometry ” subsection at the beginning of the “ Methods ” section (page 8, fist paragraph).
Biological cells play an important role in our modeling approach, since our approach is not strictly continuum based but combines continuous finite strain models with discrete cellular approaches. To clarify our approach, we have changed the following points correspondingly in order to make it clearer:
We have added an additional Figure (Fig. 6 a ) to the manuscript where we graphically present the relationship of biological versus numerical cells;

we have added a subsection to the manuscript (“ Combining continuous and discrete tissue models ”), where we provide detailed information (page 8);

we have added additional information to the “Active deformations” subsection (page 11); and

we have added a corresponding sentence to the abstract.

 2.
The model (like all mathematical models) chooses specific functional forms for, e.g., diffusing factor production and active force response to the diffusing factor. What results do the authors think are general, and robust to changes in the specific form of Eq. 3 or Eq. 4? What specific qualitative (or quantitative) predictions are made? For example, do the authors propose that, in general, straindependent feedback does not give gastrulationlike invaginations, while stressmediated feedback does? Do the authors propose that, in general, thinner tissues will have more bumps, as in Fig. 2? In the absence of general conclusions, it is difficult to evaluate the impact of this work.
We are a little bit cautious with general conclusions, since we know that patterns resulting from simulation of tissues restricted to the 2D plane may distinctly differ from those of full 3D tissues, for several reasons (examples are give with the second paragraph, page 5). However, we have added information regarding the robustness ob observed results at several places within the manuscript, e.g. page 4 paragraph 2 (robustness regarding initial conditions), page 4 paragraph 3 (robustness regarding the type of active deformations), page 5, last sentences (robustness regarding the mechanical invariant used within the feedback), and page 5, first paragraph (nonrobustness of pattern symmetry). Finally, we have added a corresponding sentence to the Conclusions section (page 7).
 3.
It is difficult to intuitively understand the three mechanical invariants (strain, compression and stress). The authors should add heat maps showing the three invariants somewhere. For example, the twelve configurations in Fig. 1 could be accompanied by duplicates with heat maps showing the invariant used. It would also greatly enhance the accessibility of the paper if the authors could add a conceptual description of the difference between strain, compression and stress in the text.
In order to address this point, we have added the following information to the manuscript:  4.
Tissue growth plays an important role in morphogenesis, especially in embryogenesis. The authors should state in the text that they are neglecting tissue growth.
We have added a corresponding statement to the manuscript (page 4, paragraph 3).
 5.
It is unusual that the authors consider inertia in their model, because biological processes at this scale are in the lowReynolds, noninertial regime. The authors should state that the results do not depend on massdensity or the inertia terms, or note if they do. (Is there a separation of timescales between extremely fast mechanics and slow diffusion, production and degradation of the diffusing factor?)
This question is related to remark 3 of referee 2 (Konstantin Doubrovinski). Inertial effects do not play a role in this model. However, instead of completely removing this term and solving for the stationary limit, we add the time derivative for reasons of stability. Temporal scales however are separated, which is realized by the parameter ε in Eq. ( 3 ). Numerically, the specific choice of ε does not change the final pattern.
 6.
Inconsistency: The main text says Fig. 2 is compressionmediated, while the caption says it is strainmediated.
We have changed this
 7.
There are several grammatical and spelling errors. For example, on p.11 Line 16, change “I turned out that?” to “It turned out that?” There is also some language usage that is unusually colloquially, notably: “Anyway, the scale of the morphogen is not crucial, since?”, suggest change to “The scale of the morphogen does not influence results, since?”
We have changed this
Revision 2
Reviewer 1: Marek Kimmel – Rice University, Houston, Texas
No additional comments on the manuscript.
Reviewer 2: Konstantin Doubrovinski (nominated by Ned Wingreen) – Universität des Saarlandes, Saarbrücken, Germany
 1.
I believe in the text following reactiondiffusion equations diffusion constants initially denoted with d’s are (erroneously?) referred to by μ’s.
We have changed this.
Reviewer 3: Jun Allard (nominated by William Hlavacek) – University of California, Irvine
 1.
The manuscript is much improved and much clearer in its description of the model. p.4 I suggest not using the term "deformation" to include tissue growth. Therefore, I suggest changing "However, other possible active deformations" to "However, other possible active processes". We have changed this.
 2.
From this version, it appears that the authors are making the assumption that active deformation only occurs at the boundaries of biological cells. This assumption is reasonable but could be more clearly stated? in its current writing, it?s easy to miss. I suggest adding a sentence on p.8 saying explicitly that active forces are only applied at the boundaries of biological cells, and again near Eq. 1 (on p.9).
Actually this is not the case, we apologize that our description was misleading. The active deformation is applied to the whole tissue body rather than only at the cell boundaries. However, the biological cell boundaries play indeed a very special role during this process, since the direction of deformation shows a jump at these boundaries. This is a result of the cytoskeleton pulling from both directions at the boundary region, the latter separating biological cells. To clarify this, we have added/changed corresponding sentences at page 8 and page 9, as suggested by the referee, and hope that this is clearer now.
 3.
Minor corrections: pg 2, second sentence of background  “pattering” change to “patterning” pg 3, last sentence  “bophysical” change to “biophysical” pg 8, second to last sentence  “cytosceleton” change to “cytoskeleton”
We have corrected this.
Abbreviations
 2D:

two dimensional
 3D:

three dimensional
 FEM:

finite element method
Declarations
Acknowledgements
M.M. and A.M.C. have been supported by the Emmy Noether Program of the Deutsche Forschungsgemeinschaft. F.B. and T.R. have been supported by the Landesstiftung BadenWürttemberg Juniorprofessorenprogramm. All authors acknowledge financial support by Deutsche Forschungsgemeinschaft and RuprechtKarlsUniversität Heidelberg within the funding programme Open Access Publishing.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Allena R, Munoz JJ, Aubry D. Diffusionreaction model for Drosophila, embryo development. Comput Methods Biomech Biomed Engin. 2013; 16(3):235–48.View ArticlePubMedGoogle Scholar
 Amar M, Goriely A. Growth and instability in elastic tissues. J Mech Phys Solids. 2005; 53:2284–319.View ArticleGoogle Scholar
 Becker R, Braack M, Dunne T, Meidner D, Richter T, Vexler B. Gascoigne 3D a finite element toolbox. 2005. http://www.numerik.unikiel.de/~mabr/gascoigne/. Accessed April 2016.
 Braybrook SA, Peaucelle A. Mechanochemical aspects of organ formation in Arabidopsis thaliana: the relationship between auxin and pectin. PLoS ONE. 2013; 8(3):e57813.View ArticlePubMedPubMed CentralGoogle Scholar
 Brouzes E, Farge E. Interplay of mechanical deformation and patterned gene expression in developing embryos. Curr Opin Genet Dev. 2004; 14:367–74.View ArticlePubMedGoogle Scholar
 Chanet S, Martin AC. Mechanical force sensing in tissues. Prog Mol Biol Transl Sci. 2014; 126:317–52.View ArticlePubMedPubMed CentralGoogle Scholar
 Cohen M, Georgiou M, Stevenson NL, Miodownik M, Baum B. Dynamic filopodia transmit intermittent deltanotch signaling to drive pattern refinement during lateral inhibition. Dev Cell. 2010; 19(1):78–89.View ArticlePubMedGoogle Scholar
 Conte V, Munoz JJ, Miodownik M. A 3D finite element model of ventral furrow invagination in the Drosophila melanogaster embryo. J Mech Behav Biomed Mater. 2008; 1:188–98.View ArticlePubMedGoogle Scholar
 Day SJ, Lawrence PA. Measuring dimensions: the regulation of size and shape. Development. 2000; 127(14):2977–87.PubMedGoogle Scholar
 Dekanty A, Milan M. The interplay between morphogens and tissue growth. EMBO Rep. 2011; 12(10):1003–10.View ArticlePubMedPubMed CentralGoogle Scholar
 Desprat N, Supatto W, Pouille PA, Beaurepaire E, Farge E. Tissue deformation modulates twist expression to determine anterior midgut differentiation in Drosophila embryos. Dev Cell. 2008; 15:470–7.View ArticlePubMedGoogle Scholar
 Farge E. Mechanotransduction in development. Curr Top Dev Biol. 2011; 95:243–65.View ArticlePubMedGoogle Scholar
 Gierer A, Berking S, Bode H, David C, Flick K, Hansmann G, Schaller H, Trenkner E. Regeneration of Hydra from reaggregated cells. Nat New Biol. 1972; 239(91):98–101.View ArticlePubMedGoogle Scholar
 Gierer A, Meinhardt H. A theory of biological pattern formation. Kybernetik. 1972; 12(1):30–9.View ArticlePubMedGoogle Scholar
 Gilbert S. Developmental Biology.Basingstoke, Hampshire, England: Palgrave Macmillan; 2013.Google Scholar
 Hiscock TW, Megason SG. Mathematically guided approaches to distinguish models of periodic patterning. Development. 2015; 142(3):409–19.View ArticlePubMedPubMed CentralGoogle Scholar
 Holzapfel GA. Nonlinear Solid Mechanics  A continuum approach for engineering.West Sussex, England: John Wiley & Sons, LTD; 2010.Google Scholar
 Howard J, Grill SW, Bois JS. Turing’s next steps: the mechanochemical basis of morphogenesis. Nat Rev Mol Cell Biol. 2011; 12(6):392–8.View ArticlePubMedGoogle Scholar
 Iskratsch T, Wolfenson H, Sheetz MP. Appreciating force and shape  the rise of mechanotransduction in cell biology. Nat Rev Mol Cell Biol. 2014; 15:825–33.View ArticlePubMedGoogle Scholar
 Jaeger J, Sharp DH, Reinitz J. Known maternal gradients are not sufficient for the establishment of gap domains in Drosophila melanogaster. Mech Dev. 2007; 124(2):108–28.View ArticlePubMedGoogle Scholar
 Junttila S, Saarela U, Halt K, Manninen A, Pärssinen H, Lecca MR, Brändli AW, SimsLucas S, Skovorodkin I, Vainio SJ. Functional genetic targeting of embryonic kidney progenitor cells ex vivo. J Am Soc Nephrol. 2015; 26(5):1126–37. doi:10.1681/ASN.2013060584. Epub 2014 Sep 8.View ArticlePubMedPubMed CentralGoogle Scholar
 Jurchenko C, Salaita KS. Lighting up the force: Investigating mechanisms of mechanotransduction using fluorescent tension probes. Mol Cell Biol. 2015; 35(15):2570–82.View ArticlePubMedPubMed CentralGoogle Scholar
 Keller R. Developmental biology. Physical biology returns to morphogenesis. Science. 2012; 338(6104):201–3.View ArticlePubMedGoogle Scholar
 Klumpers DD, Zhao X, Mooney DJ, Smit TH. Cell mediated contraction in 3d cellmatrix constructs leads to spatially regulated osteogenic differentiation. Integr Biol (Camb). 2013; 5(9):1174–83.View ArticleGoogle Scholar
 Koch A, Meinhardt H. Biological pattern fomation: from basic mechanisms to complex structures. Rev Mod Phys. 1994; 66:1481–507.View ArticleGoogle Scholar
 Kondo S, Miura T. Reactiondiffusion model as a framework for understanding biological pattern formation. Science. 2010; 329(5999):1616–20.View ArticlePubMedGoogle Scholar
 Kramerov AA, Golub AG, Bdzhola VG, Yarmoluk SM, Ahmed K, Bretner M, Ljubimov AV. Treatment of cultured human astrocytes and vascular endothelial cells with protein kinase ck2 inhibitors induces early changes in cell shape and cytoskeleton. Mol Cell Biochem. 2011; 349(12):125–37.View ArticlePubMedPubMed CentralGoogle Scholar
 Li B, Li F, Puskar KM, Wang JHC. Spatial patterning of cell proliferation and differentiation depends on mechanical stress magnitude. J Biomech. 2009; 42:1622–7.View ArticlePubMedPubMed CentralGoogle Scholar
 Mammoto A, Huang S, Ingber DE. Filamin links cell shape and cytoskeletal structure to rho regulation by controlling accumulation of p190RhoGAP in lipid rafts. J Cell Sci. 2007; 120(Pt 3):456–67.View ArticlePubMedGoogle Scholar
 Mammoto A, Mammoto T, Ingber DE. Mechanosensitive mechanisms in transcriptional regulation. J Cell Sci. 2012; 125(Pt 13):3061–73.View ArticlePubMedPubMed CentralGoogle Scholar
 Mammoto T, Ingber DE. Mechanical control of tissue and organ development. Development. 2010; 137(9):1407–20.View ArticlePubMedPubMed CentralGoogle Scholar
 Mammoto T, Mammoto A, Ingber DE. Mechanobiology and developmental control. Annu Rev Cell Dev Biol. 2013; 29:27–61.View ArticlePubMedGoogle Scholar
 MarciniakCzochra A. ReactiondiffusionODE models of pattern formation. In: Evolutionary Equations with Applications in Natural Sciences, Lecture Notes in Mathematics. Heidelberg, Germany: Springer: 2015.Google Scholar
 Martin AC, Goldstein B. Apical constriction: themes and variations on a cellular mechanism driving morphogenesis. Development. 2014; 141(10):1987–98.View ArticlePubMedPubMed CentralGoogle Scholar
 Mendez MG, Kojima SI, Goldman RD. Vimentin induces changes in cell shape, motility, and adhesion during the epithelial to mesenchymal transition. FASEB J. 2010; 24(6):1838–51.View ArticlePubMedPubMed CentralGoogle Scholar
 Mercker M, Hartmann D, MarciniakCzochra A. A mechanochemical model for embryonic pattern formation: coupling tissue mechanics and morphogen expression. PLoS ONE. 2013; 8(12):e82617.View ArticlePubMedPubMed CentralGoogle Scholar
 Mercker M, Köthe A, MarciniakCzochra A. Mechanochemical symmetry breaking in Hydra aggregates. Biophys J. 2015; 108(9):2396–407.View ArticlePubMedPubMed CentralGoogle Scholar
 Munoz JJ, Barrett K, Miodownik M. A deformation gradient decomposition method for the analysis of the mechanics of morphogenesis. J Biomech. 2007; 40(6):1372–80.View ArticlePubMedGoogle Scholar
 Murray J. Mathematical Biology II: Spatial models and biomedical applications. Heidelberg: SpringerVerlag Berlin; 2003.Google Scholar
 Nakamasu A, Takahashi G, Kanbe A, Kondo S. Interactions between zebrafish pigment cells responsible for the generation of turing patterns. Proc Natl Acad Sci U S A. 2009; 106(21):8429–34.View ArticlePubMedPubMed CentralGoogle Scholar
 Nelson CM. Geometric control of tissue morphogenesis. Biochim Biophys Acta. 2009; 1793(5):903–910.View ArticlePubMedPubMed CentralGoogle Scholar
 Nelson CM, Jean RP, Tan JL, Liu WF, Sniadecki NJ, Spector AA, Chen CS. Emergent patterns of growth controlled by multicellular form and mechanics. Proc Natl Acad Sci. 2005; 102:11594–9.View ArticlePubMedPubMed CentralGoogle Scholar
 Okuda S, Inoue Y, Watanabe T, Adachi T. Coupling intercellular molecular signalling with multicellular deformation for simulating threedimensional tissue morphogenesis. Interface Focus. 2015; 5(2):20140095.View ArticlePubMedPubMed CentralGoogle Scholar
 Patwari P, Lee RT. Mechanical control of tissue morphogenesis. Circ Res. 2008; 103:234–43.View ArticlePubMedPubMed CentralGoogle Scholar
 Philipp I, Aufschnaiter R, Ozbek S, Pontasch S, Jenewein M, Watanabe H, Rentzsch F, Holstein TW, Hobmayer B. Wnt/betacatenin and noncanonical wnt signaling interact in tissue evagination in the simple eumetazoan Hydra. Proc Natl Acad Sci U S A. 2009; 106(11):4290–5.View ArticlePubMedPubMed CentralGoogle Scholar
 Prelich G. Gene overexpression: uses, mechanisms, and interpretation. Genetics. 2012; 190(3):841–54.View ArticlePubMedPubMed CentralGoogle Scholar
 Rodriguez EK, Hoger A, McCulloch AD. Stressdependent finite growth in soft elastic tissues. J Biomech. 1994; 27(4):455–67.View ArticlePubMedGoogle Scholar
 Rodriguez M, McGarry P, Sniadecki N. Review on cell mechanics: Experimental and modeling approaches. Appl Mech Rev. 2013; 65:060801–35.View ArticleGoogle Scholar
 Rogers KW, Schier AF. Morphogen gradients: from generation to interpretation. Annu Rev Cell Dev Biol. 2011; 27:377–407.View ArticlePubMedGoogle Scholar
 Sawyer JM, Harrell JR, Shemer G, SullivanBrown J, RohJohnson M, Goldstein B. Apical constriction: a cell shape change that can drive morphogenesis. Dev Biol. 2010; 341(1):5–19.View ArticlePubMedPubMed CentralGoogle Scholar
 Shyer AE, Tallinen T, Nerurkar NL, Wei Z, Gil ES, Kaplan DL, Tabin CJ, Mahadevan L. Villification: how the gut gets its villi. Science. 2013; 342(6155):212–8.View ArticlePubMedPubMed CentralGoogle Scholar
 Soriano J, Rüdiger S, Pullarkat P, Ott A. Mechanogenetic coupling of Hydra symmetry breaking and driven turing instability model. Biophys J. 2009; 96(4):1649–60.View ArticlePubMedPubMed CentralGoogle Scholar
 Turing AM. The chemical basis of morphogenesis. Phil Trans R Soc London B. 1953; 237:37–72.View ArticleGoogle Scholar
 Turner N, Grose R. Fibroblast growth factor signalling: from development to cancer. Nat Rev Cancer. 2010; 10(2):116–29.View ArticlePubMedGoogle Scholar
 Urdy S. On the evolution of morphogenetic models: mechanochemical interactions and an integrated view of cell differentiation, growth, pattern formation and morphogenesis. Biol Rev Camb Philos Soc. 2012; 87(4):786–803.View ArticlePubMedGoogle Scholar
 Vignaud T, Blanchoin L, Thery M. Directed cytoskeleton selforganization. Trends Cell Biol. 2012; 22(12):671–82.View ArticlePubMedGoogle Scholar
 Wozniak MA, Chen CS. Mechanotransduction in development: a growing role for contractility. Nat Rev Mol Cell Biol. 2009; 10(1):34–43.View ArticlePubMedPubMed CentralGoogle Scholar
 Wyczalkowski MA, Chen Z, Filas BA, Varner VD, Taber LA. Computational models for mechanics of morphogenesis. Birth Defects Res C Embryo Today. 2012; 96(2):132–52.View ArticlePubMedPubMed CentralGoogle Scholar