Pattern formation by two-layer Turing system with complementarysynthesis

Journal of Theoretical Biology 322 (2013) 33–45

Contents lists available at SciVerse ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Pattern formation by two-layer Turing system with complementary synthesis Hironori Fujita a,n, Masayoshi Kawaguchi a,b a b

Division of Symbiotic Systems, National Institute for Basic Biology, National Institute for Natural Sciences, Okazaki 444-8585, Japan Department of Basic Biology, School of Life Science, Graduate University for Advanced Studies (SOKENDAI), Okazaki 444-8585, Japan

H I G H L I G H T S c c c c c

We study a two-layer Turing system with complementary synthesis of molecules The Turing condition of this system is determined by a linear stability analysis The system needs strong regulatory interactions of molecules for the Turing condition Pattern formation by the system is examined in fixed or expanding 2D space We use this system to explain pattern formation in the shoot meristem of plants

a r t i c l e i n f o

abstract

Article history: Received 9 August 2012 Received in revised form 6 December 2012 Accepted 11 January 2013 Available online 20 January 2013

Many multicellular organisms have a layered structure. The interaction between these layers plays an essential role in many developmental processes, and key molecules involved in these processes are often expressed in a layer-specific manner. On the other hand, pattern formation of organisms has been frequently discussed in connection with the Turing system. However, the Turing system has so far been studied mainly in single-layered space. In this paper, we thus investigate a two-layer Turing system with complementary synthesis, in which two interacting molecules are exclusively synthesized in different layers. From a linear stability analysis, we determine the Turing condition of the complementary system, and show that this condition requires stronger regulatory interactions of the molecules than that of the system with usual ubiquitous synthesis. We then confirm that this complementary system affects pattern types in fixed and expanding two-dimensional spaces in a similar way to the system with ubiquitous synthesis. In addition, the two-layer system includes two types of diffusion, lateral and transversal, and these have distinct effects on pattern formation with lateral diffusion mainly determining the periodicity of patterns generated and transversal diffusion affecting pattern type. These results suggest that the transversal diffusion functions as a time delay in the two-layer system. Finally, we apply this complementary system to explain pattern formation of the shoot apical meristem of plants. These findings provide an understanding of pattern formation caused by the interaction between cell layers in multicellular organisms. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Reaction–diffusion system Turing pattern Two-layer system Complementary synthesis Shoot apical meristem

1. Introduction Living organisms frequently show self-organized pattern formation during their developmental processes. This issue has been discussed in connection with the reaction–diffusion system that was first proposed by Turing (1952) (Meinhardt, 1982; Murray, 2003; Kondo and Miura, 2010). This system consists of diffusible molecules interacting with each other, and can generate stable patterns, namely Turing patterns, even from a homogeneous field.

n

Corresponding author. Tel.: þ81 564 55 7563; fax: þ 81 564 55 7564. E-mail address: [email protected] (H. Fujita).

0022-5193/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jtbi.2013.01.008

On the other hand, many multicellular organisms have a layered structure, and it is well known that interaction among these layers plays a crucial role in many processes of cell differentiation and organ development. In vertebrates, the interaction between the overlying ectoderm and the underlying mesenchymal cells is required for the normal development of many organs such as the limb bud (Zeller et al., 2009), lung ¨ 2006), kidney (Costantini and Kopan, 2010), (Cardoso and Lu, gustatory papilla (Jung et al., 2004), tooth (Salazar-Ciudad and Jernvall, 2010), mammary gland, and hair follicle (Mikkola and Millar, 2006; Sick et al., 2006; Schneider et al., 2009). In addition, key molecules involved in these developmental events are often expressed in a layer-specific manner. For example, the initial

34

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

stage of hair follicle formation requires the interaction of Wnt and its inhibitor Dkk, which are specifically expressed in the overlying epithelial layer and the underlying mesenchymal cells, respectively (Sick et al., 2006). Therefore, to understand morphogenesis of multicellular organisms, it is important to understand pattern formation caused by the interaction between layers. Pattern formation by the Turing system has so far been studied mainly in single-layered space. Conversely, studies on layered structure have developed in the last decade (Yang et al., 2002; Yang and Epstein, 2003, 2004; Ji and Li, 2005, 2006; Kytta¨ et al., 2007; Vasquez et al., 2008; Liu et al., 2010; Catlla´ et al., 2012). These studies showed that coupled layer systems can lead to unusual patterns not observed in single-layered systems, such as superposition patterns combining to form stripes and/or spots through the interaction of two Turing modes with different wavelengths (Yang et al., 2002; Yang and Epstein, 2003). In addition, patterns caused by coupling layers have also been experimentally investigated in chemical systems such as the chlorine dioxide–iodine–malonic acid reaction and the Belousov–Zhabotinsky reaction (Berenstein et al., 2004; Epstein et al., 2008; Mı´guez et al., 2011). However, most of these researches concentrate on the case where each layer alone can generate stable patterns with its corresponding pattern mode. In contrast, few theoretical studies exist on complementary multiple-layer systems in which molecules are synthesized in a layer-specific manner. In such systems, pattern formation never occurs in each layer alone but is caused only after interaction among layers. One of the most interesting pattern formations caused by the interaction between cell layers is the morphogenesis of the shoot apical meristem (SAM) of plants. The SAM is located in the apex of shoots, and contains a population of undifferentiated cells that have high cell division activity and constantly produce lateral organs. Extensive molecular studies revealed that the SAM

morphogenesis is centrally governed by the feedback regulation between WUSCHEL (WUS) and CLAVATA3 (CLV3) (Fig. 1c) (Stahl and Simon, 2010; Barton, 2010; Perales and Reddy, 2012). WUS, encoding a transcriptional factor with a homeodomain, moves between cells and positively regulates the CLV3 expression by directly binding to its promoter (Mayer et al., 1998; Yadav et al., 2011). In addition, WUS also activates its own expression by a positive feedback loop (Leibfried et al., 2005; Gordon et al., 2009). Conversely, CLV3, a small secreted peptide containing a CLE domain, has an inhibitory effect on WUS expression via the transmembrane receptors of CLV1, CLV2, CORYNE, and RPK2 (Clark 1995; Fletcher and Brand (1999); Brand et al., 2000; ¨ Muller et al., 2008; Kinoshita et al., 2010). The regulatory relationship between WUS and CLV3 is very similar to that of the activator–inhibitor system, one of the most famous reaction– diffusion systems (Meinhardt, 1982). Interestingly, WUS and CLV3 show complementary synthesis in a cell layer-dependent manner, where CLV3 is expressed in the outermost cell layers of the SAM and WUS below them. Several mathematical models of the SAM regulation have been reported based on the feedback regulation of WUS and CLV ¨ (Jonsson et al., 2005; Nikolaev et al., 2007; Hohm et al., 2010). Moreover, we previously proposed a reaction–diffusion model that incorporates spatial restriction and area expansion by cell division, and showed that this model can successfully explain many experimental observations (Fujita et al., 2011). We therefore proposed that the reaction–diffusion mechanism plays an essential role in the SAM pattern formation. However, since our previous model is implemented with a single cell layer, WUS and CLV3 are synthesized in the same cells and the effect of the separation between their expressions is not incorporated into the model. Accordingly, it remains to be resolved what effect this complementary synthesis has on the SAM pattern formation.

Fig. 1. Models for two-layer systems. Schematic representation of two-layer activator–inhibitor systems with (a) usual ubiquitous synthesis and (b) complementary expression of u (activator) and v (inhibitor). (c) The SAM of plants is centrally controlled by the feedback regulation between WUS and CLV3, which are highly expressed in the CZ but are exclusively synthesized in distinct layers; CLV3 is synthesized in the outermost cell layers but WUS in the underlying cell layers (the OC). (d) Schematic representation of the two-layer SAM model in the presence of the complementary synthesis of WUS and CLV3, implemented with the spatial restriction of WUS synthesis.

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

Thus, we here investigate the complementary two-layer Turing system, in which two interacting molecules are exclusively synthesized in different layers. First, we perform a linear stability analysis to determine whether or not this system has the ability to generate stable patterns. We next examine what effect this complementary synthesis has on pattern formation in fixed or expanding two-dimensional space. Finally, we confirm whether this complementary system can explain SAM pattern formation as in our previous model.

2. Model 2.1. Two-layer reaction–diffusion model We consider the reaction–diffusion system containing two components u and v in two coupled layers of one- or twodimensional space (Fig. 1a and b). The two components regulate each other and diffuse along each layer (lateral diffusion) and between overlapping two layers (transversal diffusion). The concentration dynamics of u and v in layer L ( ¼1, 2) is described by the differential equations   @uL du 2 ¼ f L ðuL ,vL Þ þ r uL þhu uL uL @t g

ð1aÞ

  @vL dv 2 ¼ g L ðuL ,vL Þ þ r vL þ hv vL vL @t g

ð1bÞ

where L indicates the layer other than L. The right-hand side of each equation contains three terms for reaction kinetics, lateral diffusion, and transversal diffusion. The coefficients of u and v are du ¼ d and dv ¼ dd in the lateral diffusion, and hu ¼h and hv ¼ Zh in the transversal diffusion. Thus, d is the ratio of v to u in the lateral diffusion, and Z is that in the transversal diffusion. In addition, g indicates a parameter with respect to spatial scale for examining the effect of area growth in Section 3.2.2. That is, after an initial calculation with g ¼ g0 for 1000.0 time units, g is increased according to

g ¼ g0 ert

ð2Þ

where g0 and r are positive constants. For the reaction dynamics, we consider the activator–inhibitor system, which is one of the most famous reaction–diffusion systems and is frequently applied to the pattern formation of living systems. We use a simple form of reaction terms with u (activator) and v (inhibitor) described by f L ðu,vÞ ¼ jðF þAs uBvÞEu and g L ðu,vÞ ¼ CuDv

ð3Þ

where As, B, C, D, E, and F are positive constants, and j is a constraint function such that j(x)¼0 if xo0, j(x)¼ jmax if x 4 jmax, and j(x) ¼x if 0 rxr jmax (Fig. 1a). The upper limit of j is set to be jmax ¼E umax ¼E u0/K, where umax is the upper limit of u, u0 is the steady state of u, and K ¼u0/umax o1 is a positive constant. This constraint condition of u synthesis results in that of u concentration: 0 ru rumax ¼u0/K. The first and second terms on the right side of each equation indicate synthesis and degradation, respectively. We then introduce the condition of complementary synthesis so that u and v are exclusively synthesized in layers L¼2 and L¼1, respectively (Fig. 1b), which is written as f L ðu,vÞ ¼ ½jðF þ As uBvÞL2 Eu and g L ðu,vÞ ¼ ½CuL1 Dv

ð4Þ

where [x]c are constraint functions that have the value x in layer L¼1 (if c ¼L1) or in layer L¼ 2 (if c¼L2), and otherwise have the value 0. Thus, this system never forms patterns in each layer alone, but can only generate patterns through the interaction between the two layers.

35

The initial values of variables are given as their steady state with random fluctuation of 1.0%. For numerical simulations, Eqs. (1a, 1b) is discretized with time step Dt, and calculations are performed with the periodic boundary condition with domain size 200.0 (grid: 200) (Figs. 2 and 4) or 200.0  200.0 (grid: 200  200) (Figs. 6 and 7). Parameter values used are described in the figure legends. 2.2. Two-layer SAM model We consider a SAM model using two overlapping cell-layers, which extends our previous model with single cell layer (Fujita et al., 2011). Numerical simulations are carried out by iterating the following four steps in order: reaction–diffusion dynamics, cell removal, cell division, and vertex dynamics. 2.2.1. Reaction–diffusion step The SAM consists of the CZ (central zone) at the center and its surrounding PZ (peripheral zone), and the outer region of the SAM called the OZ (organ zone) (Fig. 1c) (Clark, 1997). The OC (organizing center) is defined as the lower region of the CZ, where WUS is highly expressed. Because the SAM of plants usually maintains its constant size despite active cell division, it is possible that the SAM is constrained by a spatial restriction to prevent abnormal proliferation of stem cells. It is thought that the CZ induces the PZ differentiation from the OZ, but the molecular mechanism is still unclear (Reddy and Meyerowitz, 2005). In our previous SAM model, therefore, we introduced a presumptive diffusible molecule z, which is produced in the CZ and has the ability to induce PZ differentiation. We here extend our previous model to include coupling two cell layers, in which CLV3 (the inhibitor, denoted by v) and WUS (the activator, denoted by u) are exclusively synthesized in the upper layer (L¼1) and lower layer (L¼2), respectively (Fig. 1d). The spatial restriction on the SAM is implemented by limiting WUS synthesis, and is therefore applied to the lower layer (L¼2). Cells of the lower layer are divided into three regions, the OC, PZ, and OZ, from the center to the periphery. The OC is defined as the cells that have high WUS concentrations (u 4uc ¼2u0). In the OC, z is specifically synthesized and diffuses to generate a gradient, inducing PZ differentiation if z is larger than the threshold zs. WUS is produced in the SAM (the OC and PZ) of the lower layer (L¼2) only. As a result, the concentrations of u, v, and z in layer L cell i obey the following equations: X duL,i ¼ ½jðG þ As uL,i BvL,i ÞSAM EuL,i þ du ðuL,j uL,i Þ dt j X þ hu ðwL,k-L,i uL,k wL,i-L,k uL,i Þ

ð5aÞ

k

X dvL,i ¼ ½CuL,i L1 DvL,i þ dv ðvL,j vL,i Þ dt j X ðwL,k-L,i vL,k wL,i-L,k vL,i Þ þ hv

ð5bÞ

k

X X dzL,i ¼ ½Z OC zL,i þdz ðzL,j zL,i Þ þ hz ðwL,k-L,i zL,k wL,i-L,k zL,i Þ dt j k ð5cÞ where Z, dz, and hz are constants. The four terms on the right-hand side of each equation represent synthesis, decay, lateral diffusion, and transversal diffusion, respectively. [x]c are constraint functions that have the value x in the SAM (z4zc) of layer L¼2 (if c¼SAM), in layer L¼1 (if c¼L1), or in the OC (u 4uc) of layer L¼2 (if c¼OC), and have the value 0 otherwise. The transversal movement of molecules across the layers from layer L cell i to

36

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

layer L cell k depends on the ratio of their overlapping area (SL,i\L,k ) to that of the exporting cell (SL,i): wL,i-L,k ¼ SL,i\L,k =SL,i . The initial values of variables are given as their steady state values with random fluctuation of 1.0%, and the reflective boundary condition is used. In every reaction–diffusion step, calculations are carried out for a total time of 50.0 and a time step of Dt.

for a total time of 100.0 and a time step of 0.1, which leads to an almost equilibrium state.

3. Results 3.1. Linear stability analysis

2.2.2. Cell removal step The cells that are considerably apart from the SAM with zozd ¼0.01 are removed in the lower layer (L¼2), while the cells that do not overlap with any lower cell are in the upper layer (L¼1). 2.2.3. Cell division step The largest cell divides perpendicular to its long axis in each of the upper and lower cell layers. Concentrations of the molecules in a divided cell are set to be the same as those of its mother cell. 2.2.4. Vertex dynamics step Vertex dynamics is carried out by a method similar to that in Nagai and Honda (2001). Each cell layer consists of polygonal cells that contact tightly with each other without gaps. Thus, each cell is defined by its surrounding vertices which follow the equation of motion, dr i =dt ¼ ri U ¼ ð@U=@xi ,@U=@yi Þ, where ri ¼(xi, yi) is the 2D-position vector of vertex i, and ri is the nabla differential operator with respect to ri. The right-hand side of the equation represents a potential force (driving force), that is, minus the gradient of the potential U. Differentiation of the potential U with respect to time P gives the inequality dU=dt ¼  ðdr i =dtÞ2 r0, indicating that veri tices move to decrease U. The potential U contains two terms related to cell edge length (Lk) and cell area (Si), which are expressed by P P vertex coordinates, U ¼ s Lk þ k ðSi Sstd Þ2 , where s ¼0.02, i k k ¼0.1, and Sstd ¼1.0 are constants. Thus, vertices move so that cell edge length becomes smaller while maintaining a constant area of each cell. In every vertex dynamics step, calculations are carried out

To examine whether or not the two-layer system with complementary synthesis can generate stable patterns by selforganization (namely Turing patterns), we perform a detailed linear stability analysis. The condition for generating Turing patterns (referred to as the Turing condition) consists of two requirements: the steady state is stable under small perturbations in the absence of lateral diffusion (Req. A), but it becomes unstable if lateral diffusion is introduced when a standing wave develops (Req. B). This type of destabilization of the steady state is known as the Turing bifurcation. 3.1.1. Two-layer system with ubiquitous synthesis Let us first consider the simple case without the expression separation, in which u and v are synthesized in both layers, described by Eqs. (1a, 1b) and (3) (Fig. 1a). The Turing condition for this case is consequently equivalent to that of the usual single-layer system, and may be summarized as Ao D, B 4 AD=C, A 4 D=d, and B oðdA þ DÞ2 =4dC

ð6Þ

where A ¼As–E and d ¼dv/du (Appendix A). When this Turing condition is satisfied, the Turing instability is caused for wave numbers in the range: k ok o k þ

ð7Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 2 2 where ¼g kc ðBCADÞ=dd Þ and kc ¼ ðdADÞ=2dd. This result is consistent with the numerical analysis, in which Turing 2 k7

2 ðkc 7

Fig. 2. Parameter conditions for generating Turing patterns in one-dimensional two-layer systems with (a) ubiquitous synthesis and (b) complementary synthesis of u and v. Circles indicate generation of Turing patterns, in which the steady state is stable without lateral diffusion and stable patterns are generated by introducing the diffusion. Triangles indicate that stable patterns are generated but the steady state is unstable without lateral diffusion. Crosses indicate no stable patterns generated. Gray regions indicate the parameter area that satisfies the Turing condition predicted by Eq. (6)(a) or by Eq. (12)(b). Calculations are performed using Eqs. (1a, 1b) and (3) in (a) or Eqs. (1a, 1b) and (4) in (b), with C¼D ¼ E¼ F¼ d ¼h ¼ g ¼1.0, d ¼ Z ¼20.0, K¼0.1, and Dt¼ 0.001.

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

patterns are observed within the parameter condition predicted by Eq. (6) (Fig. 2a). 3.1.2. Two-layer system with complementary synthesis We next examine the influence of the complementary synthesis given by Eqs. (1a, 1b) and (4), in which the activator (u) is specifically synthesized in layer L¼ 2 and the inhibitor (v) in layer L¼1 (Fig. 1b). We consider a small perturbation around the steady state, uL ¼ uL0 þ U L eltikx and vL ¼ vL0 þ V L eltikx , where UL and VL are perturbation amplitudes, and l and k are the growth rate with time and the wave number of the perturbation, respectively. This perturbation is introduced to the linearized equations of Eq. (1a, 1b) around the steady state, and we then 2 obtain the following relationship: ðJlI þ k DL =gÞw ¼ O, where w ¼(U1, U2, V1, V2)T 0 B B J¼B B @

Ehu

hu

0

hu

Ahu

0

C

0

Dhv

0

0

hv

0

1

0

du

B C B 0 C C, and DL ¼ B B 0 C @ A Dhv 0

0

0

B

du

0

hv

0

dv

0

0

0

1

C 0 C C 0 C A dv

To have a non-zero solution for w, the determinant of the matrix must be zero to obtain the dispersion relation: 2

detðlÞ ¼ 9JlI þ k DL=g 9 ¼ 0

ð8Þ

where det(l) is a quartic function of l. 3.1.2.1. In the absence of lateral diffusion (k ¼0). The steady state is stable if Re(l) is negative for all solutions of Eq. (8). In the absence of lateral diffusion (or k¼ 0), this requirement (i.e. Req. A) is consequently summarized as B 4 B0 ¼

DðD þ 2ZhÞððE þhÞAEhÞ

ð9aÞ

Zh2 C 2

B o Bp ¼

2ð2hþ EAÞðZh þDÞðha þ bÞð4Zð1 þ ZÞh þ ð1 þ2ZÞha þ bÞ

Zh2 ða þ 2ð1 þ ZÞhÞ2 C ð9bÞ

37

where a ¼ 2D þEA and b ¼ ðD þ EÞðDAÞ (for details see Appendix B). The inequalities in Eq. (9) can be expressed in other forms with h4h0 and h4hp, respectively, where h0 is the largest solution of B ¼B0 and hp is the largest solution of B ¼Bp. Thus, Eq. (9) (i.e. Req. A) can be also described in the form: h4 maxðh0 ,hp Þ

ð10Þ

Since h represents the transversal diffusivity, this inequality indicates that stronger interactions between the layers are required for the steady state to be stable in the absence of lateral diffusion.

3.1.2.2. In the presence of lateral diffusion (ka0). We next consider the condition causing Turing instability (i.e. Req. B), that is, the steady state being destabilized by introducing lateral diffusion causing a standing wave is developed. This condition is satisfied if one of the solutions of Eq. (8) is positive for some k40. This is equivalent to det(0)¼a0(k2, h) being negative for some k40, where 2 2 0 2 0 2 0 4 0 2 a0 ðk ,hÞ ¼ BC Zh þðdd k þDÞðdd k þ D þ 2ZhÞðd k þd ð2hA0 Þk  0 0 0 A hEAÞ, A ¼A E, and d ¼d/g. Therefore, this condition depends on the sign of a0(k2, h), and is then determined by the zero contour line of a0(k2, h)¼0, which is referred to here as h¼ha0(k2) (Fig. 3). That is, the steady state is stable if h4ha0(k2) corresponding to a0(k2, h)40, but becomes unstable if hoha0(k2) corresponding to a0(k2, h)o0. This contour line runs through (k2, h)¼(0, h0) and (A/d0 , 0), and is classified into the following three cases:

(i) ha0(k2) decreases monotonically as k2 increases (Fig. 3(i)); (ii) ha0(k2) has a maximum point (denoted by k¼km and h¼hm) (Fig. 3(ii)), which is determined by the equations:  @a0  2 a0 ðkm ,hm Þ ¼ 0 and ¼0 ð11Þ  2 @ðk Þk2 ¼ k2 ,h ¼ hm m

(iii) ha0(k2) increases for kok  and decreases for k4k þ as k2 increases, and a0(k2, h)o0 for any h4 0 in k  okok þ

Fig. 3. The Turing condition in the two-layer system with complementary synthesis described by Eqs. (1a, 1b) and (4). The dotted area indicates the parameter condition that never satisfies Eq. (9) (i.e. Req. A) for any h 40. The Turing condition depends on the zero contour line of a0(k2, h) ¼ 0, namely h ¼ ha0(k2). The shape of this contour line is divided into cases (i), (ii), and (iii). The gray area in the A-B plane indicates the Turing condition predicted by Eq. (12). A1 is the A value of the intersection of h ¼h0 and h ¼ hp, and Am is that of h ¼h0 and h ¼ hm (see Appendices B and D).

38

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

(Fig. 3(iii)), where k þ and k  are described below in the next section. Therefore, with respect to a0(k2, h), the condition causing the Turing instability needs a0(k2, h)4 0 at k¼0 and a0(k2, h) o0 for some k 40. This condition is satisfied when hm 4h 4h0 in case (ii) or h4h0 in case (iii), whereas Turing instability never happens in case (i). Taking into account Eq. (10) (i.e. Req. A), the Turing condition is then defined as hm 4h4max(h0, hp) in case (ii) or h4max(h0, hp) in case (iii), and is consequently represented by the parameter area enclosed by the three lines of h¼h0 (corresponding to Eq. (9a)), h¼hp (Eq. (9b)), and h¼hm (Eq. (11)) in Fig. 3. As a result, the Turing condition may also be summarized as hm 4h 4 maxðho ,hp Þ

ð12Þ

This result is confirmed by numerical simulations, in which Turing patterns are observed in the parameter area predicted by Eq. (12) (Fig. 2b). These results indicate that the two-layer system with complementary synthesis has the ability to generate Turing patterns similar to the single layer. However, the Turing condition for this complementary system (Fig. 2b, gray area) is shifted towards larger A and B compared with that of the usual ubiquitous synthesis (Fig. 2a, gray area). Because A and B indicate regulatory strengths of self-enhancement of the activator and negative feedback between the activator and inhibitor, respectively, this result suggests that the complementary system requires stronger regulatory interaction between the molecules for generating Turing patterns than the ubiquitous synthesis system.

Fig. 4. Effect of transversal diffusivity on pattern formation in one-dimensional two-layer system with complementary synthesis. Circles and crosses are the same as described in the legend of Fig. 2. The gray area indicates the Turing condition predicted by Eq. (12). Calculations are performed using Eqs. (1a, 1b) and (4) with A¼ 1.2, B¼3.0, C¼D ¼ E¼ F¼ d ¼ g ¼ 1.0, d ¼ 20.0, K¼ 0.1, and Dt ¼0.001.

3.1.2.3. In the infinite limit of transversal diffusivity. Let us consider the Turing condition in the infinite limit of h, which occurs in case (iii) of the previous section. This condition is consequently summarized as A0 o D0 , B 4A0 D0 =C, A0 4D0 =d, and B o ðdA0 þD0 Þ2 =4dC

ð13Þ

0

where A ¼ AE and D0 ¼2D (for details see Appendix C). If this requirement is satisfied, Turing instability is caused in the wave qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 4 2 number k  okok þ , where k 7 ¼ gðkc 7 kc ðBCA0 D0 Þ=dd Þ and 2

kc ¼ ðdA0 D0 Þ=2dd. The Turing condition given by Eq. (13) coincides with that of the usual activator–inhibitor system, described by Eq. (6), except for a transformation of parameters. This result suggests that this complementary two-layer system behaves like a single-layer system under the condition of large h, that is, strong diffusivity across the layers. However, the parameter transformation shifts the Turing condition towards increasing A and B, indicating that the complementary system needs stronger regulatory interaction between molecules to cause Turing instability than the usual system with ubiquitous synthesis. 3.1.2.4. Effects of lateral and transversal diffusivities. The lateral and transversal diffusivities are represented by d ¼du and h ¼hu, respectively, with d ¼dv/du being the ratio of the lateral diffusion of the inhibitor to the activator and Z ¼hv/hu being that of the transversal diffusion. The two diffusions affect the Turing condition in different manners. For example, it is known that the Turing condition requires d 41 (i.e. dv 4du) in the usual single layer system (Murray, 2003), and this inequality is also true in the complementary two-layer system (Appendix D). In contrast, the transversal diffusivity h also affects the Turing condition, but a restriction on its ratio Z is not necessarily required (Fig. 4). That is, Turing patterns can be generated even if Z o1 (i.e. hv o hu).

Fig. 5. Effect of lateral and transversal diffusivities on the wave number with the fastest growth rate (km). The value of km is calculated by numerically solving Eq. (11) with A ¼1.7, B ¼3.0, C¼ D ¼E ¼F ¼ g ¼ 1.0, d ¼ Z ¼20.0, and K¼ 0.1. The Turing condition is satisfied when h 4max(h0, hp), where h0 E2.17 and hp E1.82.

In addition, the two diffusions have distinct effects on the spatial periodicity of the patterns generated. That is, the wave number with the largest growth rate (km) decreases (that is, the spatial period increases) as the lateral diffusivity (d) increases, while the transversal diffusivity (h) does not have a large effect on the wave number (Fig. 5).

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

3.2. Pattern formation in two-dimensional space Since the complementary two-layer system has the ability to form stable patterns autonomously, we investigate its effect on pattern types generated in fixed or expanding two-dimensional areas.

3.2.1. Pattern types generated in fixed area It is known that stable patterns generated by the reaction– diffusion dynamics in 2D space are divided into three types: spot patterns (Fig. 6a), stripe patterns (Fig. 6b), and reverse spot patterns (Fig. 6c) (Shoji et al., 2003). In the activator–inhibitor system, it is also known that these pattern types depend on the relative distance between the upper and lower limits of the steady state of the activator (Shoji et al., 2003). Thus, in the case of the two-layer system, we examine the effect of this relative position (K ¼u0/umax) in layer L¼2 where u is synthesized. In both ubiquitous and complementary synthesis, pattern types generated change drastically in the order of the spot

39

pattern, stripe pattern, and reverse spot pattern as K increases, that is, as the steady state becomes closer to the upper limit (Fig. 6d and e). Conversely, the pattern type is not affected much by B (a parameter of the reaction dynamics) (Fig. 6d and e) or by h (transversal diffusivity) (Fig. 6f). These results suggest that the pattern type is mainly affected by the constraint condition of the activator, but not by the transversal diffusivity (h).

3.2.2. Effect of area expansion It has often been reported that generated patterns respond during growth in one-dimensional (Crampin et al., 1999, 2002; Murray, 2003; Gaffney and Monk, 2006; Neville et al., 2006; Seirin-Lee and Gaffney, 2010; Seirin-Lee et al., 2011) and two-dimensional space (Meinhardt, 1982; Varea et al., 1997; Madzvamuse et al., 2003; Plaza et al., 2004; Gjorgjieva and Jacobsen, 2007). Furthermore, extensive studies have been made of models with area growth and shape change used to explain experimental observations such as patterns in sea shells (Meinhardt, 2003), fish skin pigmentation (Kondo and Asai,

Fig. 6. Pattern types generated in fixed two-dimensional space are divided into three groups: (a) spot patterns (filled circles); (b) stripe patterns (squares); and (c) reverse spot patterns (open circles). Ubiquitous synthesis is implemented in (d) using Eqs. (1a, 1b) and (3), and complementary synthesis in (a–c, e, f) using Eqs. (1a, 1b) and (4). (a–c) Darker gray indicates higher levels of the activator in layer L¼ 2. (d–f) Crosses indicate that no stable patterns are generated. The parameters are set to be A ¼ 0.8 (d) or 1.0 (a–c, e, f), B ¼4.0 (a–c) or 3.5 p (f), C ¼D ¼ E¼ F¼ d ¼1.0, d ¼ Z ¼20.0, h ¼ 1.0 (a–e), K¼ 0.2 (a), 0.35 (b), or 0.45 (c), Dt ¼0.003125 (d) or 0.0025 (a–c, e, f), and ffiffiffiffiffiffiffiffiffiffiffi g ¼ k20 =k2m0 , where k2m0 ¼ ððA þ DÞ þ ðd þ 1Þ BC=dÞ=ðd1Þd, k0 ¼ 2np=200:0, and n¼ 10.0 (d) or 15.0 (a–c, e, f). (f) The Turing condition is satisfied when hm 4 h 4max(h0, hp), where h0 E 0.595, hp E0.456, and hm E1.38.

40

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

1995, Venkataraman et al., 2011), vertebrate limb skeletons (Hentschel and Glimm (2004); Miura et al., 2006; Zhu et al., 2010), teeth (Salazar-Ciudad and Jernvall, 2004, 2010), green algae (Harrison et al., 1981, 2001), and plant shoot apexes (Holloway and Harrison, 2008; Fujita et al., 2011). However, the effects of area growth in twodimensional space are less studied than those in one dimension and in particular, parameter dependency on pattern formation has not been examined in detail. We thus examine the effect of area expansion with respect to the spot pattern generated in a fixed two-dimensional area. In the cases of both ubiquitous and complementary synthesis, the manner in which the spot pattern changes during area expansion is classified into the following four proliferation modes.

(i). Elongation mode: each initial spot continues to elongate during area expansion (Fig. 7d).

(ii). Division mode: each spot initially elongates and is then equally divided to give rise to two parts, resulting in generating multiple spots (Fig. 7c). (iii). Emergence mode: spots move apart from each other by area expansion, and new spots emerge from between existing spots, leading to multiple spots as in the division mode (Fig. 7b). (iv). Fluctuation mode: this mode generates multiple spots with relatively weak activator concentrations. During area expansion, these spots behave in a fluctuating manner where spots move, are divided into two spots, merge with neighboring spots, or fade away (Fig. 7a). In contrast to the pattern formation in fixed area that is affected mainly by K, the proliferation mode depends on both K and B. However, these parameters have different effects. For example, as B increases, the proliferation type shifts in the order

Fig. 7. Pattern proliferation caused by area expansion is divided into four types: (a) fluctuation mode (diamonds); (b) emergence mode (triangles); (c) division mode (squares); and (d) elongation mode (circles). Ubiquitous synthesis is implemented in (e) using Eqs. (1a, 1b)–(3), and complementary synthesis in (a–d, f, g) using Eqs. (1a, 1b), (2) and (4). (a–d) Darker gray indicates higher levels of the activator in layer L ¼2. (e–g) Crosses indicate that no stable patterns are generated. The parameters are set to be A¼ 0.8 (e) or 1.0 (a–d, f, g), B¼ 5.4 (a), 4.4 (b), 3.4 (c), 2.6 (d), or 3.5 (g), C¼ D¼ E ¼F ¼d ¼ 1.0, d ¼ Z ¼ 20.0, h ¼1.0 (a–f), K¼ 0.1 (a-d), Dt ¼0.000625 (e) or 0.0008 pffiffiffiffiffiffiffiffiffiffiffi 2 2 2 (a–d, f, g), r ¼0.001, and g ¼ k0 =km0 , where km0 ¼ ððA þ DÞ þ ðd þ 1Þ BC=dÞ=ðd1Þd, k0 ¼ 2np=200:0, and n¼ 4.5 (e) or 6.0 (a–d, f, g). (g) The Turing condition is satisfied when hm 4h 4max(h0, hp), where h0 E0.595, hp E0.456, and hm E1.38.

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

41

of elongation mode, division mode, emergence mode, and fluctuation mode (Fig. 7e and f). Conversely, K affects the pattern selection between the division and emergence modes with the division mode tending to occur as K increases, but only a slight effect on parameter conditions of the elongation and fluctuation modes (Fig. 7e and f). That is, the elongation and fluctuation modes are strongly affected by B but only slightly by K, whereas the division and emergence modes are affected by both parameters. In addition, the transversal diffusivity (h¼hu) affects the proliferation mode in a similar way to B (Fig. 7g), possibly because h exerts its effect on the proliferation type by changing the reaction dynamics. These results suggest that the complementary two-layer system can generate stable patterns as seen in the usual single layer Turing system. 3.3. SAM pattern formation The SAM of plants is a typical example of complementary synthesis. That is, in A. thaliana, CLV3 (inhibitor) is produced in the outer cell layers, but the expression of WUS (activator) is limited to the underlying cell layers (Fig. 1c). However, it remains unclear what effect this expression separation has on SAM pattern formation. On the other hand, we have previously developed a SAM model that successfully explains many experimental results, but this model is implemented with a single cell layer (Fujita et al., 2011). Thus, to clarify the effect of complementary synthesis between WUS and CLV3 on SAM pattern formation, we extend our previous model to include complementary synthesis using a simplified two-layered cell network (Fig. 1d), and then examine whether this extended model can explain SAM pattern formation or not. The extended model developed here described by Eq. (5) can generate various SAM patterns as in the case of our previous model (Fujita et al., 2011). That is, SAM patterns generated by the model are classified into the following six groups depending on the proliferation mode of the OC and the strength of the spatial restriction of the SAM (Fig. 8).

(i) Homeostasis pattern: a single OC spot continues to be located at the center of the SAM with a constant SAM size (Fig. 9c). This pattern corresponds to the wild-type SAM of A. thaliana. (ii) Fasciation pattern: an initial OC spot continues to extend in the elongation mode, resulting in the SAM containing an extremely elongated OC (Fig. 9a). This morphological feature is similar to the fasciation that is observed in the field environment and in many mutants such as clv mutants (Iliev and Kitin, 2011; Leyser and Furner, 1992; Clark, 1993). (iii) Fluctuation pattern: OC spots with weak WUS expression multiply in the fluctuation mode, leading to a SAM with many weak OC spots (Fig. 9d). This pattern is related to the phenotype of the wus mutant (Laux et al., 1996), because the SAM of wus shows a patchy pattern of the WUS promoter expression with very weak levels compared to those of the wild type (Fujita et al., 2011). (iv) Dichotomous pattern: an OC spot is elongated and divided into two parts in division mode followed by division of the PZ, resulting in two independent SAMs (Fig. 9b). This is associated with the dichotomous branching of shoots, which is observed in plants such as mistletoe (Viscum album) and in mutants such as clv mutants (Leyser and Furner, 1992) and klv of Lotus japonicus (Oka-Kira et al., 2005, Miyazawa et al., 2010). (v) Monopodial pattern: a new OC spot gives rise to the PZ in the emergence mode, and then two SAMs are formed by splitting

Fig. 8. SAM pattern formation in two-cell-layer system with complementary expression of WUS and CLV3. SAM patterns generated by the model are classified into six groups: fasciation pattern (filled circles); dichotomous pattern (open squares); homeostasis pattern (open circles); fluctuation pattern (diamonds), monopodial pattern, and multiplication pattern in the emergence mode (filled triangles) or in the division mode (filled squares). Calculations are performed using Eq. (5), and the parameters are set to be A ¼ 0.6, C¼ D ¼E ¼ F¼ dz ¼hz ¼1.0, d¼ h ¼ 0.25, d ¼ Z ¼ 20.0, K¼ 0.1, zc ¼ 0.1, zd ¼ 0.02, and Dt ¼0.025.

of the PZ as in the case of the dichotomous pattern. This pattern morphologically resembles the monopodial branching of shoots, but no corresponding examples are detected in A. thaliana. It should be noted that the monopodial branching (or lateral branching) observed in many plants is thought to be caused by the feedback regulation between auxin and its ¨ transmembrane carrier PIN1 (Smith et al., 2006; Jonsson et al., 2006). (vi) Multiplication pattern: OC spots continue to multiply in the division or emergence mode as a result of a chain reaction between the OC proliferation and SAM expansion, giving rise to a massive SAM containing multiple OC spots. This morphological feature is similar to that of the amp1/pt mutant (Chaudhury et al., 1993; Vidaurre et al., 2007). These SAM patterns are affected by B (corresponding to the CLV3 function) and Z (associated with the spatial restriction) (Fig. 8). For example, under the condition of moderately tight spatial restriction (for instance, Z¼2.0 in Fig. 8), as B increases the SAM pattern shifts in the order: fasciation pattern, dichotomous pattern, homeostasis pattern, and fluctuation pattern (Fig. 9). This pattern change predicted by the model is consistent with the experimental results where clv mutants (corresponding to B reduction) show fasciated or bifurcated stems (Leyser and Furner, 1992; Clark 1993) and overexpression of CLV3 gene (corresponding to B increasing) causes a wus-like phenotype (Lenhard and Laux, 2003). On the other hand, the homeostasis pattern changes to the multiplication pattern by increasing Z, that is, by reducing the spatial restriction strength (Fig. 8). These results are similar to

42

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

CLV3

WUS

B = 6.0

B = 8.0

B = 10.0

B = 14.0

Fig. 9. Examples of time evolution of the SAM in (a) fasciation pattern, (b) dichotomous pattern, (c) homeostasis pattern, and (d) fluctuation pattern. WUS in the lower cell layer (L ¼2) is indicated in red, and CLV3 in the upper cell layer (L¼ 1) is in blue. The OZ region of the lower layer is indicated in green. Calculations are performed using Eq. (5), and the parameters used are the same as those of Fig. 8 except for Z¼ 3.0 and B¼6.0 (a), 8.0 (b), 10.0 (c), or 14.0 (d).

those of our previous model using a single layer, and are consistent with many experimental observations (for a detailed discussion see Fujita et al., 2011). In conclusion, these results show that the complementary two-layer system can explain SAM pattern formation.

4. Discussion Many multicellular organisms have a layered structure, and communications between layers are known to be very important for various pattern formations. In the case of vertebrates, many superficial organs essentially involve the interaction between an epithelium layer and the underlying mesenchymal tissue (Jung et al., 2004; Mikkola and Millar, 2006; Schneider et al., 2009). In addition, essential molecules involved in these developmental processes are often expressed in a layer-specific manner. Therefore, understanding the pattern formation caused by layerspecific synthesis is very important for understanding the development of multicellular organisms. However, few theoretical studies have examined such complementary multiple-layer systems that molecules are synthesized in a layer-specific manner. Thus, in this paper, we have investigated the complementary two-layer system, and showed that this system can generate pattern formation similar to the single layer system. This finding provides an understanding of pattern formation caused by the interaction between cell layers in multicellular organisms. In biological systems, molecules are usually synthesized through gene transcription and mRNA translation processes that consist of many kinetic steps, and hence these steps can be influential as a time delay in regulatory dynamics. It is well known that time delays have various effects on pattern formation in Turing systems (Li and Ji, 2004; Ji and Li, 2005; Gaffney and Monk, 2006; Hu et al., 2008; Sen et al., 2009; Seirin-Lee and Gaffney, 2010; Seirin-Lee et al., 2010). For example, a time delay can induce or suppress Turing patterns depending on its strength, inducing Turing instability when the undelayed system remains homogeneously stable, while a further increase in delay causes temporal oscillations to Turing patterns (Fig. 10a) (Li and Ji, 2004; Ji and Li, 2005; Sen et al., 2009; Seirin-Lee et al., 2010). Because molecules in our models are assumed to be synthesized by

transcription-translation reactions, it is expected that the effects of time delay described above are also true for our two-layer systems. On the other hand, the two-layer system with complementary synthesis described by Eqs. (1a, 1b) and (4) is equivalent to a single-layer system with four diffusible molecules (u1, u2, v1, and v2), in which u2 and v1 correspond to the activator and inhibitor, respectively, and signal transduction from u2 to v1 is mediated by u1 and from v1 to u2 by v2 (Fig. 10b). It is thus expected that these signaling steps function as time delay, and the strength of which depends on the transversal diffusivity (h¼ hu), because this diffusivity affects the regulatory strength between u1 and u2 and between v1 and v2. In this report, we showed that the transversal diffusion can affect the Turing pattern formation, such that Turing patterns are generated for an intermediate strength of transversal diffusion (hm 4h4max(h0, hp)), but turn to a spatially homogeneous state if the diffusivity is strong (h 4hm), while weak diffusivity (homax(h0, hp)) causes divergence or oscillatory behavior (Fig. 10a). This effect of transversal diffusivity is opposite to that caused by time delays as described in the above. This result is consistent because the effect of time delays becomes stronger as the transversal diffusivity decreases. As a result, this evidence suggests that the complementary synthesis between cell layers has an effect as a time delay in multiple-layer systems. One of the most interesting pattern formations caused by the interaction between layers is the SAM of plants. Our previous model explains the SAM patterning well by using the reaction– diffusion mechanism based on the regulatory relationship between WUS and CLV3 (Fujita et al., 2011). We therefore proposed that the reaction–diffusion mechanism plays a central role in the SAM pattern formation. However, it is not known what influence the separate expression of WUS and CLV3 has on SAM pattern formation, since the previous model is implemented with a single cell layer. The study in this paper shows that this expression restriction does not interfere with SAM pattern formation. Therefore, this result supports our argument that the reaction–diffusion mechanism is indispensable for SAM pattern formation. Our model showed that the SAM pattern can be formed regardless of the expression separation of WUS and CLV3. This raises the question of why WUS and CLV3 are expressed in different cell layers but not in the same. The reason may be

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

43

vL, leading to the reduced equations: @u0 du 2 @v0 dv 2 ¼ f u u0 þf v v0 þ r u0 and ¼ g u u0 þg v v0 þ r v0 @t g @t g

ðA:1Þ

where u0 ¼ ðu1 þu2 Þ=2u0 and v0 ¼ ðv1 þv2 Þ=2v0 , and fu, fv, gu, and gv are the partial differentials evaluated at the steady state. Because Eq. (A.1) is independent of h and Z, pattern formation occurs regardless of transversal diffusion. Eq. (A.1) represents a reaction–diffusion system with two components, and the Turing condition for this case is well known: fu þgv o0, fugv  fvgu 40, dvfu þdugv 40, and (dvfu þdugv)2  4dudv(fugv fvgu)40 (Murray, 2003). In the case of the reaction terms of Eq. (3), this Turing condition may be summarized as Eq. (6).

Appendix B. Turing condition in the absence of lateral diffusion In the absence of lateral diffusion (or k¼0), Eq. (8) becomes 4

Fig. 10. Effects of transversal diffusivity and time delay in Turing system. (a) A time delay can induce or suppress Turing patterns depending on its strength. That is, it induces Turing instability when the undelayed system remains homogeneously stable, and further increases to the delay cause temporal oscillations to Turing patterns. Conversely, the transversal diffusivity in the two-layer system has an opposite effect to the time delay, because the effect of the time delay becomes stronger as the transversal diffusivity decreases. (b) The two-layer system with complementary synthesis described by Eqs. (1a, 1b) and (4) is equivalent to the single-layer system where u2 and v1 are the activator and inhibitor, respectively, and signal transduction between the two molecules is mediated by u1 and v2. These signaling processes can function as a time delay, where its strength depends on the transversal diffusivity.

because of the adaptive strategy of plants rather than developmental constraints. That is, the SAM is a vital organ for plants, because it contains stem cells producing all the tissues for aerial parts of plants including reproductive organs. In addition, it is likely that the most important region for maintaining the SAM is the OC where WUS is highly expressed, because the SAM may be regenerated even after physical damage if the OC maintains its function whereas a defect in the OC may lead directly to the SAM breakdown. Thus, it is very important for plants to protect the OC from the external environment. Therefore, a plant with the OC located in an inner region of the SAM may be more adaptive than one with the OC in an outer region. Plants may have evolved such that the OC is located in inner SAM. In conclusion, this study showed that the two-layer system with complementary synthesis can generate Turing patterns, but that it requires stronger regulatory interaction between the molecules to cause Turing instability than does the usual ubiquitous system. In addition, this system affects pattern formation in fixed and expanding two-dimensional space in a similar way to the case with usual ubiquitous synthesis. Using this complementary two-layer system, the SAM pattern formation is successfully explained.

Acknowledgments This research was supported by Grants-in-Aid for Scientific Research (22128006) from the Ministry of Education, Culture, Sports, Science and Technology of Japan.

3

2

detðlÞ ¼ a4 l þ a3 l þ a2 l þ a1 l þ a0 ¼ 0

ðB:1Þ

where 2

a0 ¼ BC Zh DðD þ2ZhÞðA0 h þEAÞ 2

a1 ¼ 2ð2DA0 ÞZh 2A0 Dð1þ ZÞh þ2D2 hA0 D2 2EAðD þ ZhÞ 2

a2 ¼ 4Zh þ2Dð2 þ ZÞhA0 ð2D þh þ 2ZhÞ þ D2 EA a3 ¼ 2DA0 þ2ð1 þ ZÞh a4 ¼ 1 and A0 ¼ AE. For the steady state to be stable, Re(l) must be negative for all the solutions of Eq. (B.1). Using the Routh– Hurwitz stability criterion, this condition can be expressed as a0 4 0,a1 40,a3 4 0, and a3 a2 a1 4a4 a1 2 þ a0 a3 2

ðB:2Þ

These inequalities may be rewritten respectively as 2

B 4 B0 , A oA1 ¼

2ð2D þ EÞZh þ 2DðD þE þ EZÞh þED2 2

2Zh þ 2ðD þ DZ þEZÞh þ DðD þ2EÞ

,

A oA3 ¼ 2D þ Eþ 2ð1 þ ZÞh, and B o Bp

ðB:3Þ

where B0 and Bp are described in Eq. (9). Since we always have A3 4A1, we can remove the condition AoA3 from Eq. (B.3). Furthermore, the three lines B¼ B0 (i.e. h¼h0), B¼Bp (i.e. h ¼hp), and A¼A1 intersect at a point (see Fig. 3). Thus, we can also remove condition AoA1, and we consequently obtain Eq. (9) as the condition of (Req. A).

Appendix C. Turing condition in the infinite limit of transversal diffusivity From the argument described in the text, the Turing condition can be summarized as a0(0, h) 40, p(0)40, and a0(k2, h)o0 for some k40 (which correspond to h4h0, h 4hp, and h ohm, 2 respectively), where pðk Þ ¼ a3 a2 a1 ða4 a2 1 þa0 a3 2 Þ. Because we 5 0 have pð0Þ ¼ h ðp þOð1=hÞÞ and p0 ¼ 16ð1 þ ZÞZ2 ðD0 A0 Þ where A0 ¼ AE and D0 ¼ 2D, the sign of p(0) coincides with that of p0 in the infinite limit of h. Thus the condition p(0) 40, that is p0 40, is satisfied if A0 o D0

Appendix A. Turing condition for two-layer system with ubiquitous synthesis After the linearization of Eq. (1a, 1b) around its steady state (u0, v0), the resulting equations are summed with respect to uL or

ðC:1Þ 2

2

2

On the other hand, since we have a0 ðk ,hÞ ¼ Zh ða00 ðk Þ þ 2 02 4 0 2 Oð1=hÞÞ and a00 ðk Þ ¼ 4dd k þ 2d ðD0 dA0 Þk þ BC A0 D0 , the sign 2 0 2 of a0(k , h) becomes the same as that of a0 (k ) in the infinite limit of h. Accordingly, the condition with respect to a0(k2, h) is equivalent to a00 (0) 40 and a00 (k2)o0 for some k40.

44

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

Consequently, this condition results in 0

0

0

0 2

B 4 A D =C, A 4 D =d, and B o ðdA þ D Þ =4dC 0

0

ðC:2Þ

As a result, the Turing condition in the infinite limit of h can be described as Eqs. (C.1) and (C.2).

Appendix D. Constraint condition on lateral diffusivity The lines h¼h0 and h¼hp intersect on the line A¼A1, and the lines h¼h0 and h¼hm intersect on A ¼ Am ¼

2

2ð2D þ EdÞZh þ 2DðD þ EZ þ EdÞh þ ED2 2dZh2 þ 2ðDd þ DZ þ EdZÞh þ DðD þ 2EdÞ

(see Fig. 3). Therefore, the parameter area satisfying the Turing condition is restricted to Am oAoA1. Thus, the existence of this parameter area requires Am oA1, which is satisfied if and only if d 41 (i.e. dv 4du). References Barton, M.K., 2010. Twenty years on: the inner workings of the shoot apical meristem, a developmental dynamo. Dev. Biol. 341, 95–113. Berenstein, I., Dolnik, M., Yang, L., Zhabotinsky, A.M., Epstein, I.R., 2004. Turing pattern formation in a two-layer system: superposition and superlattice patterns. Phys. Rev. E 70, 046219. Brand, U., Fletcher, J.C., Hobe, M., Meyerowitz, E.M., Simon, R., 2000. Dependence of stem cell fate in Arabidopsis on a feedback loop regulated by CLV3 activity. Science 289, 617–619. ¨ J., 2006. Regulation of early lung morphogenesis: questions, Cardoso, W.V., Lu, facts and controversies. Development 133, 1611–1624. Catlla´, A.J., McNamara, A., Topaz, C.M., 2012. Instabilities and patterns in coupled reaction-diffusion layers. Phys. Rev. E 85, 026215. Chaudhury, A.M., Letham, S., Craig, S., Dennis, E.S., 1993. amp1 - a mutant with high cytokinin levels and altered embryonic pattern, faster vegetative growth, constitutive photomorphogenesis and precocious flowering. Plant J. 4, 907–916. Clark, S.E., 1997. Organ formation at the vegetative shoot meristem. Plant Cell 9, 1067–1076. Clark, S.E., Running, M.P., Meyerowitz, E.M., 1993. CLAVATA1, a regulator of meristem and flower development in Arabidopsis. Development 119, 397–418. Clark, S.E., Running, M.P., Meyerowitz, E.M., 1995. CLAVATA3 is a specific regulator of shoot and floral meristem development affecting the same processes as CLAVATA1. Development 121, 2057–2067. Costantini, F., Kopan, R., 2010. Patterning a complex organ: branching morphogenesis and nephron segmentation in kidney development. Dev. Cell 18, 698–712. Crampin, E.J., Gaffney, E.A., Maini, P.K., 1999. Reaction and diffusion on growing domains: scenarios for robust pattern formation. Bull. Math. Biol. 61, 1093–1120. Crampin, E.J., Gaffney, E.A., Maini, P.K., 2002. Mode-doubling and tripling in reaction-diffusion patterns on growing domains: a piecewise linear model. J. Math. Biol. 44, 107–128. Epstein, I.R., Berenstein, I.B., Dolnik, M., Vanag, V.K., Yang, L., Zhabotinsky, A.M., 2008. Coupled and forced patterns in reaction-diffusion systems. Philos. Transact. A Math. Phys. Eng. Sci. 366, 397–408. Fletcher, J.C., Brand, U., Running, M.P., Simon, R., Meyerowitz, E.M., 1999. Signaling of cell fate decisions by CLAVATA3 in Arabidopsis shoot meristems. Science 283, 1911–1914. Fujita, H., Toyokura, K., Okada, K., Kawaguchi, M., 2011. Reaction-diffusion pattern in shoot apical meristem of plants. PLoS One 6, e18243. Gaffney, E.A., Monk, N.A., 2006. Gene expression time delays and Turing pattern formation systems. Bull. Math. Biol. 68, 99–130. Gjorgjieva, J., Jacobsen, J., 2007. Turing patterns on growing spheres: the exponential case. Discrete Contin. Dyn. Syst. Ser. B, 436–445. Gordon, S.P., Chickarmane, V.S., Ohno, C., Meyerowitz, E.M., 2009. Multiple feedback loops through cytokinin signaling control stem cell number within the Arabidopsis shoot meristem. Proc. Natl. Acad. Sci. USA 106, 16529–16534. Harrison, L.G., Wehner, S., Holloway, D.M., 2001. Complex morphogenesis of surfaces: theory and experiment on coupling of reaction-diffusion patterning to growth. Faraday Discuss. 120, 277–294. Harrison, L.G., Snell, J., Verdi, R., Vogt, D.E., Zeiss, G.D., Green, B.R., 1981. Hair morphogenesis in Acetabularia mediterranea: temperature-dependent spacing and models of morphogen waves. Protoplasma 106, 211–221. Hentschel, H.G.E., Glimm, T., Glazier, J.A., Newman, S.A., 2004. Dynamical mechanisms for skeletal pattern formation in the vertebrate limb. Proc. R. Soc. B 271, 1713–1722. Hohm, T., Zitzler, E., Simon, R., 2010. A dynamic model for stem cell homeostasis and patterning in Arabidopsis meristems. PLoS One 5, e9189. Holloway, D.M., Harrison, L.G., 2008. Pattern selection in plants: coupling chemical dynamics to surface growth in three dimensions. Ann. Bot. 101, 361–374.

Hu, H.X., Li, Q.S., Ji, L., 2008. Superlattice patterns and spatial instability induced by delay feedback. Phys. Chem. Chem. Phys. 10, 438–441. Iliev, I., Kitin, P., 2011. Origin, morphology, and anatomy of fasciation in plants cultured in vivo and in vitro. Plant Growth Regul. 63, 115–129. Ji, L., Li, Q.S., 2005. Turing pattern formation in coupled reaction-diffusion system with distributed delays. J. Chem. Phys. 123, 94509. Ji, L., Li, Q.S., 2006. Turing pattern formation in coupled reaction-diffusion systems: effects of sub-environment and external influence. Chem. Phys. Lett. 424, 432–436. Jung, H.S., Akita, K., Kim, J.Y., 2004. Spacing patterns on tongue surface-gustatory papilla. Int. J. Dev. Biol. 48, 157–161. ¨ Jonsson, H., Heisler, M.G., Shapiro, B.E., Meyerowitz, E.M., Mjolsness, E., 2006. An auxin-driven polarized transport model for phyllotaxis. Proc. Natl. Acad. Sci. USA 103, 1633–1638. ¨ Jonsson, H., Heisler, M., Reddy, G.V., Agrawal, V., Gor, V., Shapiro, B.E., Mjolsness, E., Meyerowitz, E.M., 2005. Modeling the organization of the WUSCHEL expression domain in the shoot apical meristem. Bioinformatics 21 (1), i232–i240. Kinoshita, A., Betsuyaku, S., Osakabe, Y., Mizuno, S., Nagawa, S., Stahl, Y., Simon, R., Yamaguchi-Shinozaki, K., Fukuda, H., Sawa, S., 2010. RPK2 is an essential receptor-like kinase that transmits the CLV3 signal in Arabidopsis. Development 137, 3911–3920. Kondo, S., Asai, R., 1995. A reaction-diffusion wave on the skin of the marine angelfish Pomacanthus. Nature 376, 765–768. Kondo, S., Miura, T., 2010. Reaction-diffusion model as a framework for understanding biological pattern formation. Science 329, 1616–1620. ¨ K., Kaski, K., Barrio, R.A., 2007. Complex turing patterns in non-linearly Kytta, coupled systems. Physica A 385, 105–114. ¨ Laux, T., Mayer, K.F., Berger, J., Jurgens, G., 1996. The WUSCHEL gene is required for shoot and floral meristem integrity in Arabidopsis. Development 122, 87–96. Leibfried, A., To, J.P., Busch, W., Stehling, S., Kehle, A., Demar, M., Kieber, J.J., Lohmann, J.U., 2005. WUSCHEL controls meristem function by direct regulation of cytokinin-inducible response regulators. Nature 438, 1172–1175. Lenhard, M., Laux, T., 2003. Stem cell homeostasis in the Arabidopsis shoot meristem is regulated by intercellular movement of CLAVATA3 and its sequestration by CLAVATA1. Development 130, 3163–3173. Leyser, H.M.O., Furner, I.J., 1992. Characterisation of three shoot apical meristem mutants of Arabidopsis thaliana. Development 116, 397–403. Li, Q.S., Ji, L., 2004. Control of Turing pattern formation by delayed feedback. Phys. Rev. E 69, 046205. Liu, F.C., He, Y.F., Pan, Y.Y., 2010. Superlattice patterns in coupled Turing systems. Commun. Theor. Phys. 53, 971–976. Madzvamuse, A., Wathen, A.J., Maini, P.K., 2003. A moving grid finite element method applied to a model biological pattern generator. J. Comput. Phys. 190, 478–500. ¨ Mayer, K.F., Schoof, H., Haecker, A., Lenhard, M., Jurgens, G., Laux, T., 1998. Role of WUSCHEL in regulating stem cell fate in the Arabidopsis shoot meristem. Cell 95, 805–815. Meinhardt, H., 1982. Models of Biological Pattern Formation. Academic press, London. Meinhardt, H., 2003. The Algorithmic Beauty of Sea Shells, Third ed. Springer, New York. ˜ uzuri, A.P., 2011. Interaction of chemical Mı´guez, D.G., Dolnik, M., Epstein, I., Mun patterns in coupled layers. Phys. Rev. E 84, 046210. Mikkola, M.L., Millar, S.E., 2006. The mammary bud as a skin appendage: unique and shared aspects of development. J. Mammary Gland Biol. Neoplasia 11, 187–203. Miura, T., Shiota, K., Morriss-Kay, G., Maini, P.K., 2006. Mixed-mode pattern in Doublefoot mutant mouse limb—Turing reaction-diffusion model on a growing domain during limb development. J. Theor. Biol. 240, 562–573. Miyazawa, H., Oka-Kira, E., Sato, N., Takahashi, H., Wu, G.J., Sato, S., Hayashi, M., Betsuyaku, S., Nakazono, M., Tabata, S., Harada, K., Sawa, S., Fukuda, H., Kawaguchi, M., 2010. The receptor-like kinase KLAVIER mediates systemic regulation of nodulation and non-symbiotic shoot development in Lotus japonicus. Development 137, 4317–4325. ¨ Muller, R., Bleckmann, A., Simon, R., 2008. The receptor kinase CORYNE of Arabidopsis transmits the stem cell-limiting signal CLAVATA3 independently of CLAVATA1. Plant Cell 20, 934–946. Murray, J.D., 2003. Mathematical Biology II: Spatial Models and Biomedical Applications, third ed. Springer, Berlin. Nagai, T., Honda, H., 2001. A dynamic cell model for the formation of epithelial tissues. Philos. Mag. B 81, 699–719. Neville, A.A., Matthews, P.C., Byrne, H.M., 2006. Interactions between pattern formation and domain growth. Bull. Math. Biol. 68, 1975–2003. Nikolaev, S.V., Penenko, A.V., Lavrekha, V.V., Melsness, E.D., Kolchanov, N.A., 2007. A model study of the role of proteins CLV1, CLV2, CLV3, and WUS in regulation of the structure of the shoot apical meristem. Russ. J. Dev. Biol. 38, 383–388. Oka-Kira, E., Tateno, K., Miura, K., Haga, T., Hayashi, M., Harada, K., Sato, S., Tabata, S., Shikazono, N., Tanaka, A., Watanabe, Y., Fukuhara, I., Nagata, T., Kawaguchi, M., 2005. klavier (klv), a novel hypernodulation mutant of Lotus japonicus affected in vascular tissue organization and floral induction. Plant J. 44, 505–515. Perales, M., Reddy, G.V., 2012. Stem cell maintenance in shoot apical meristems. Curr. Opin. Plant. Biol. 15, 10–16. ˜ o, F., Padilla, P., Barrio, R.A., Maini, P.K., 2004. The Plaza, R.G., Sa´nchez-Gardun effect of growth and curvature on pattern formation. J. Dyn. Diff. Eq. 16, 1093–1121.

H. Fujita, M. Kawaguchi / Journal of Theoretical Biology 322 (2013) 33–45

Reddy, G.V., Meyerowitz, E.M., 2005. Stem-cell homeostasis and growth dynamics can be uncoupled in the Arabidopsis shoot apex. Science 310, 663–667. Salazar-Ciudad, I., Jernvall, J., 2004. How different types of pattern formation mechanisms affect the evolution of form and development. Evol. Dev. 6, 6–16. Salazar-Ciudad, I., Jernvall, J., 2010. A computational model of teeth and the developmental origins of morphological variation. Nature 464, 583–586. Schneider, M.R., Schmidt-Ullrich, R., Paus, R., 2009. The hair follicle as a dynamic miniorgan. Curr. Biol. 19, R132–R142. Seirin-Lee, S., Gaffney, E.A., 2010. Aberrant behaviours of reaction diffusion selforganisation models on growing domains in the presence of gene expression time delays. Bull. Math. Biol. 72, 2161–2179. Seirin-Lee, S., Gaffney, E.A., Monk, N.A., 2010. The influence of gene expression time delays on Gierer–Meinhardt pattern formation systems. Bull. Math. Biol. 72, 2139–2160. Seirin-Lee, S., Gaffney, E.A., Baker, R.E., 2011. The dynamics of Turing patterns for morphogen-regulated growing domains with cellular response delays. Bull. Math. Biol. 73, 2527–2551. Sen, S., Ghosh, P., Riaz, S.S., Ray, D.S., 2009. Time-delay-induced instabilities in reaction-diffusion systems. Phys. Rev. E 80, 046212. Shoji, H., Iwasa, Y., Kondo, S., 2003. Stripes, spots, or reversed spots in twodimensional Turing systems. J. Theor Biol. 224, 339–350. Sick, S., Reinker, S., Timmer, J., Schlake, T., 2006. WNT and DKK determine hair follicle spacing through a reaction–diffusion mechanism. Science 314, 1447–1450. Smith, R.S., Guyomarc’h, S., Mandel, T., Reinhardt, D., Kuhlemeier, C., Prusinkiewicz, P., 2006. A plausible model of phyllotaxis. Proc. Natl. Acad. Sci. USA 103, 1301–1306.

45

Stahl, Y., Simon, R., 2010. Plant primary meristems: shared functions and regulatory mechanisms. Curr. Opin. Plant Biol. 13, 53–58. Turing, A.M., 1952. The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. Ser. B 237, 37–72. Varea, C., Aragon, J.L., Barrio, R.A., 1997. Confined Turing patterns in growing systems. Phys. Rev. E 56, 1250–1253. Venkataraman, C., Sekimura, T., Gaffney, E.A., Maini, P.K., Madzvamuse, A., 2011. Modeling parr-mark pattern formation during the early development of Amago trout. Phys. Rev. E 84, 041923. Vidaurre, D.P., Ploense, S., Krogan, N.T., Berleth, T., 2007. AMP1 and MP antagonistically regulate embryo and meristem development in Arabidopsis. Development 134, 2561–2567. ¨ Yadav, R.K., Perales, M., Gruel, J., Girke, T., Jonsson, H., Reddy, G.V., 2011. WUSCHEL protein movement mediates stem cell homeostasis in the Arabidopsis shoot apex. Genes Dev. 25, 2025–2030. Yang, L., Epstein, I.R., 2003. Oscillatory Turing patterns in reaction–diffusion systems with two coupled layers. Phys. Rev. Lett. 90, 178303. Yang, L., Epstein, I.R., 2004. Symmetric, asymmetric, and antiphase Turing patterns in a model system with two identical coupled layers. Phys. Rev. E 69, 026211. Yang, L., Dolnik, M., Zhabotinsky, A.M., Epstein, I.R., 2002. Spatial resonances and superposition patterns in a reaction–diffusion model with interacting Turing modes. Phys. Rev. Lett. 88, 208303. Zeller, R., Lo´pez-Rı´os, J., Zuniga, A., 2009. Vertebrate limb bud development: moving towards integrative analysis of organogenesis. Nat. Rev. Genet. 10, 845–858. Zhu, J.F., Zhang, Y.T., Alber, M.S., Newman, S.A., 2010. Bare bones pattern formation: a core regulatory network in varying geometries reproduces major features of vertebrate limb development and evolution. PLoS One 5, e10892.