Gene network models robust to spatial scaling and noisy input

Mathematical Biosciences 237 (2012) 1–16

Contents lists available at SciVerse ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Gene network models robust to spatial scaling and noisy input Heather Hardway Department of Mathematics and Statistics, and Center for BioDynamics, Boston University, 111 Cummington Street, Boston, MA 02215, United States

a r t i c l e

i n f o

Article history: Received 10 August 2011 Received in revised form 12 January 2012 Accepted 1 March 2012 Available online 17 March 2012 Keywords: Reaction–diffusion Pattern formation Bicoid–Hunchback Gene network

a b s t r a c t Many biological systems are inherently noisy, yet demonstrate robustness to perturbations and changes in external influences. Such is the case in the Bicoid–Hunchback (Bcd–Hb) system, which is critical to axis specification in the developing Drosophila embryo. We use this system as motivation to explore the larger problem of how precise patterning can be achieved under imprecise conditions. While evidence suggests Bicoid gradients are uncorrelated with respect to embryo length, downstream genes, such as Hb, are expressed in a precise manner with regard to position along the anterior–posterior (AP)-axis. In addition to precision under variability of embryo length, Hb also exhibits robustness to perturbations to the regulatory network, gene dosage, and temperature. Understanding the reduced variability of patterns in this system is of interest to both experimentalists and theoreticians, lending itself well to the field of mathematical modeling. In this paper, a class of reaction–diffusion models is presented, which produce precise patterns, despite receiving noisy input and other perturbations to the system. An essential property of the network includes the existence of a strong inhibitor for the Hb representative, where the strength of the inhibition is directly related to the amount of variation that can be tolerated. With a higher inhibitory effect, larger perturbations of Bcd can be made with relatively small changes to the location of the Hb boundary. Network topology and interaction strength are the essential properties of the minimal model giving rise to the robust features, and possible interpretations are made with regard to the Bcd–Hb system. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction From what was originally one cell, a fully developed organism displays a range of cell types, each exhibiting unique characteristics. While all cells contain an identical genome, different cell types display distinct patterns of gene expression, and these differences drive cellular specification and differentiation in developing embryos as a function of position and time. For species to survive, development must be robust to both external and internal variables, such as differences in environmental conditions or gene expression levels. The fruit fly, Drosophila melanogaster, has been used as a model system in genetics and developmental biology, and its gene network has become one of the most studied and well understood. This is especially true in specification of the anterior–posterior and dorsal–ventral axis, the first steps in spatially organizing the developing embryo [1–3]. Central to the current understanding of such spatial organization during Drosophila embryonic patterning is the concept of morphogen gradients, which act to spatially distinguish cell types by forming a concentration gradient that controls the expression of other genes depending on the level of

E-mail address: [email protected] 0025-5564/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mbs.2012.03.004

morphogen at that point [4–6]. While the initial concept of morphogen gradients in distinguishing cell fates was theoretical, Drosophila soon provided the first characterized example, Bicoid (Bcd), shown to regulate the expression of downstream gap genes in a concentration-dependent fashion [7–9]. With asymmetric expression from initial mRNA localization, Bicoid forms an anterior-to-posterior concentration gradient [10]. Hunchback (Hb), positively regulated by Bcd, subsequently is expressed with high levels at the anterior and a sharp boundary near the middle of the embryo, crucial for subsequent developmental decisions and AP-axis specification [11–14]. While experimental studies have identified the hierarchal structure of the Drosophila gene network and the relationship between genes, recent studies have indicated a robust patterning of genes downstream from Bcd. Hb was shown to be highly correlated with respect to embryo length, forming a sharp border at the midpoint of the embryo, despite receiving noisy positional information form its upstream regulator, Bcd [15]. More so, the precision was also robust to perturbation of the gene network and temperature variation [15,16]. In addition to scaling properties within members of the same species, [17] showed that patterns of gene expression also scaled across species of Dipterans, with considerable variability in mean embryo length. Subsequent studies have supported the premise that Hb and other gap genes effectively filter noisy in-

2

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

put from Bcd, resulting in scaled patterns with higher precision than seen in maternal gradients [18–20]. Reduced variability of gap gene expression is also seen over time, with higher variability at early stages of development [21,22]. However, recent work has also supported the contrary, that Bcd is sufficiently precise to encode observed Hb expression [23–25]. While debate on the amount of variability in the Bcd–Hb system continues, the larger question of how precision is achieved in developmental systems remains.

cAA ¼

kA caa ; cAB ¼ cAA j; RA

cBA ¼

RA kB ccAA ; cBB ¼ jcBA : RB kA

1.1. Specific model

There exists a positively invariant region for Eq. (1), given by

A system of two interacting proteins, A and B, is considered, where each can diffuse, decay, and influence their own production and the production of the other. The embryo is modeled as a one dimensional line segment, 0 6 x 6 L, due to the homogeneity of AP patterns along the circumference of the embryo. The pairwise effect of protein i on the production of protein j is represented by the constant cij . The function f(x) represents an external gradient that is allowed to influence the production of each protein. First, the equations for A and B are rescaled to dimensionless variables, see Appendix A. The equations are further simplified by requiring the front solution to be located at the midpoint of space and robust to perturbations of f(x), see Appendix B. Therefore, we consider the following specific two-component reaction–diffusion model, viewing a and b as shifted, dimensionless protein concentrations,

R ¼ fða; bÞj  1=2 6 a 6 1=2; 1=ð2lÞ 6 b 6 1=ð2lÞg:

@a @2a ¼ Da 2  a þ gða  lcb þ qf ðxÞÞ @t @x @b @2b ¼ Db 2  lb þ gðjða  lcb þ qf ðxÞÞÞ @t @x 0 6 x 6 1; t P 0

ð1Þ

ajt¼0 ¼ a0 ðxÞ; bjt¼0 ¼ b0 ðxÞ @a @b ; j ¼0 j @x x¼0;1 @x x¼0;1 with

gðzÞ ¼

1 1  : 1 þ ecaa z 2

Eq. (1) is a specific example of models used in [26,27]. g represents a shifted, scaled production function of a linear combination of a and b. More generally, one may take g to be monotone increasing in a rough network model, where the cumulative affect of all transcription factors is assumed to be represented by the weighted sum of gene expression levels [28]. However, in this paper, we will assume g to be a sigmoidal function, i.e. g is bounded, monotone increasing, switches concavity from up to down at z ¼ 0, and is odd. In addition, we assume a constant affect of each protein on the production of another. While there is evidence to suggest some genes have dual regulatory affects and subsequent models for this [29,30], these will not be considered in this paper. The parameters are defined by

cFA kA ; cAA RA 2r k c¼1þ a A; cAA RA

q¼

Da ¼ DA =ðkA L2 Þ; k l¼ B; kA

In relation to the original, dimensional variables, A and B,

a¼

kA 1 kA 1 A ;b ¼ B : 2 2l RA RB

Provided the initial data is in R; then ðaðx; tÞ; bðx; tÞ 2 R for all x 2 ½0; 1; t P 0. a will be thought to represent a simplified Hb-like protein, b a second protein, and f a Bcd-like protein. Due to the activation of Bcd on hb and self-activation of hb, we will assume caa ; cfa > 0 [9,11,31]. See Appendix A for details regarding the original dimensional variables and parameters. The network topology indicated by Eq. (1) with c; j > 0 (Network Type I) is depicted in Fig. 1. Similarly, there are three other distinct network topologies. However, by a change of variables, network types 1 and 4 are equivalent, as are network types 2 and 3. Not surprisingly, network type 3, consisting of all positive interactions, can never produce stable, steady-state patterns. Therefore, network type 1 is the focus of the paper, and we will assume c; j > 0. The spatially homogeneous solutions of Eq. (1) in the absence of f (i.e. q ¼ 0) solve the ordinary differential equation system

v 0 ¼ v þ gðv  lcwÞ; w0 ¼ lw þ gðjðv  lcwÞÞ:

ð2Þ

Eq. (2) can exhibit monostability or bistability. The point (0, 0) is always an equilibrium, and it is the only equilibrium if ð1  cjÞg 0 ð0Þ < 1. To see this, define n :¼ v  lcw. Then steadystate solutions of Eq. (2) satisfy

n ¼ gðnÞ  cgðjnÞ: n ¼ 0 is the only solution of Eq. if ð1  cjÞg 0 ð0Þ < 1. Moreover, it is stable as a fixed point of (2), and via the classical Turing mechanism, can become unstable in Eq. (1) with unequal diffusion rates, Da – Db [32,33]. In contrast, (0, 0) is only one of three equilibria when ð1  cjÞg 0 ð0Þ > 1. (0, 0) is always unstable, while the other two are stable homogeneous states. Representative phase portraits for the kinetic equations in the monostable and bistable regimes are shown in frames (a) and (b) of Fig. 2. In contrast, the dynamics of representative solutions of the full PDE system, Eq. (1), are qualitatively similar, as illustrated in frames (c), (e) and (d), (f) of Fig. 2. Numerical simulations suggest that if Db is sufficiently large, the monotone solutions are the only stable steady states (Fig. 3). However, when the difference in diffusion rates is decreased, both monotone and non-monotone patterns persist for large times, dependent on initial conditions (Fig. 4). Numeric approximations for the solutions were obtained using a backward Euler scheme with spatial and time steps Dx ¼ 0:005; Dt ¼ 0:01. Smaller step sizes did not significantly change the approximations. 1.2. Main results and organization

Db ¼ DB =ðkA L2 Þ;

where DA ; RA ; kA ; ra are the dimensional diffusion, maximal production, decay, and basal production rates for protein A, respectively. Similarly, the dimensional interaction coefficients are defined by

In this paper, a class of reaction–diffusion models, Eq. (1), is presented that mimics the key characteristics of the Bcd–Hb system, in that the component, a, representing Hb, exhibits a sharp boundary at the midpoint of space, and this is robust to changes in (a) length of space, (b) initial conditions, and (c) functional inputs f(x), representing Bcd and other known regulators. As opposed to

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

Fig. 1. Schematic of network topology for Eq. (1) for c > 0; j > 0, i.e. Network Type 1; cij denotes the pairwise influence on the production of protein i by protein j.

previous models, the dynamics do not require a specific underlying structure to the model, i.e. bistability, a Turing instability, or a threshold mechanism. Instead, properties (a) and (b) emerge from the stability properties of the resulting shadow system, while (c) is related to the inhibitory strength in the network, which is proportional to the amount of variation that can be effectively filtered. The shadow system associated to Eq. (1), which is obtained by setting b equal to its spatial average, is analysed in Section 2. The model is shown to be robust to both changes in the length of space and initial conditions. Section 3 presents the linear analysis of the full coupled system, Eq. (1), again in the absence of f. While multiple steady-state patterns are observed, like the shadow system, the characteristic eigenmode corresponds to the monotone solutions, as long as the diffusion rates are significantly different. Section 4 presents numeric evidence that both the shadow system and full coupled system exhibit a buffering capacity to variability in external input, f(x). Of particular importance is the magnitude of the inhibitory interactions in the network. In addition, using a Heaviside approximation for g, we show that the same cannot be true for a scalar reaction–diffusion equation, modeling a single protein influenced by itself and an external gradient. This shows that at least two interacting proteins are necessary for robust pattern formation. Lastly, the results are discussed with possible implication for the Bcd–Hb system in Section 5.

1.3. Relation to literature The Bcd–Hb system has proven to be of interest across disciplines, lending itself to a theoretical framework as another means of study. Many proposed mathematical models for biological pattern formation are based on the early pioneering theories of Wolpert and Turing. In Wolpert’s French flag model, the absolute value of a morphogen determines how space is divided [4]. This model, which can also be referred to as a threshold model, is strongly dependent on the values of the morphogen at each point. If a noisy gradient is used instead, the French flag model will result in a readout with wavy, ill defined stripes [34]. Therefore, a simple threshold model cannot capture one of the key characteristics of early Bcd–Hb system in Drosophila development. However, a variation on this theme is considered in models requiring a posterior gradient, which is a reflection of the Bcd gradient, where the intersec-

3

tion of the two designates the midpoint of the embryo and location for Hb boundary [35–37]. Turing proposed a pattern formation mechanism based on two morphogens, a slow diffusing activator and fast diffusing inhibitor [38]. While contrary to intuition, the different diffusion rates can destabilize a constant equilibrium in the well-mixed system, allowing small perturbations to amplify and create spatially inhomogeneous solutions [32,33]. Linear stability analysis predicts the dominant unstable eigenmode and corresponding characteristic wavelength for pattern formation [38]. Two classic models, by Gierer-Meinhardt and Gray-Scott, have been extensively studied and shown to produce complicated spatiotemporal patterns [39–42]. Significant interest has been devoted to the scale invariance of Turing structures, in part due to the application of such models for biological phenomena, including morphogenesis. Proposed methods for fixing the characteristic wavelength under changes in total size include introducing a dependence of kinetic and diffusion coefficients on size, and similar modifications along these themes [43–45]. Similar size invariance issues were studied with regard to the growth and cutting of domains, as indicated from other developmental systems [46]. Other models were based on the hierarchical structure of the gene systems and increasing frequency of patterns in Drosophila morphogenesis [46–53]. Alternatively, spatial pattern formation can arise from multiple stable steady-state concentration levels of the morphogens. When two steady states are observed, the system is termed bistable and has been applied to various biological systems, including Drosophila patterning system [19,54–57]. Such bistable, or multistable systems, have been used to explain the sharpness of patterns and amplifications of initial discontinuities [57,58].

2. Analysis of shadow system In this section, we analyze the shadow system associated to Eq. (1) with q ¼ 0,

at ¼ Da axx  a þ gða  lcbÞ Z 1 b0 ¼ lb þ gðjða  lcbÞdx

ð3Þ

0

with boundary conditions ax ð0; tÞ; ax ð1; tÞ ¼ 0. Here, b denotes the spatial average of bðxÞ, and the prime indicates the derivative with respect to time. This shadow system is obtained from Eq. (1) by considering the regime in which b diffuses very quickly. See [59,60] for other examples of shadow systems arising in coupled reaction–diffusion equations. Time independent solutions of Eq. (3) satisfy the following differential equation coupled to an integral constraint

0 ¼ Da axx  a þ gða  lcbÞ; Z 1 0 ¼ lb þ gðjða  lcbÞdx

ð4Þ

0

with boundary conditions ax ð0Þ; ax ð1Þ ¼ 0. We will analyze both the constant and nontrivial spatially dependent solutions of Eq. (4), as well as the bifurcation structure of such solutions. We will prove the following theorem: Theorem 2.1. Let ðas ; bs Þ designate a solution to Eq. (4). 1. (Existence of Homogeneous Steady State Solutions) For Da > Dc :¼ ðg 0 ð0Þ  1Þ=p2 , only the constant steady-state solutions exist. If the system is in the monostable regime, i.e. if

4

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

0.5

0.3

0.3

0.1

0.1

w

w

0.5

−0.1

−0.1

−0.3

−0.3 −0.5

−0.5 −0.5

−0.3

−0.1

v

0.1

0.3

−0.5

0.5

−0.3

−0.1

v

(a) a(x,t)

0.5

a(x,t)

0.5

50

40

40

30

t

30 0

t 20

20

10

10

0

0.2

0.4

0.6

0.8

−0.5

1

0

0

(d) 0.5 0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0.1

0.2

0.2

0.3

0.3

0.4

0.4 0.4

0.6

0.8

1

0.8

1

0

0.1

0.2

0.6

(c) 0.4

0

0.4

x

0.5

0.5

0.2

x

a

a

0.3

(b)

50

0

0.1

0.5

0

0.2

0.4

x

x

(e)

(f)

0.6

0.8

1

Fig. 2. (a), (b) Phase plane portraits for Eq. (2); blue lines represent sample solution trajectories, maroon and orange lines the nullclines; solid circles represent stable equilibria, unfilled circles unstable equilibria; (a) caa ¼ 8; l; c; j ¼ 1, there exists a single, stable equilibrium; (b) caa ¼ 8; l; c ¼ 1; j ¼ 1=4, there exist three equilibrium, two of which are stable; (c), (d) contour plots for aðx; tÞ in Eq. (1), Da ¼ 0:01; Db ¼ 1; b0 ðxÞ ¼ 0; a0 ðxÞ ¼ 0:1  x; caa ¼ 8; q ¼ 0; (c) l; c; j ¼ 1 (as in (a)); (d) l; c ¼ 1; j ¼ 1=4 (as in (b)); (e), (f) aðx; TÞ corresponding to (c), (d) for T ¼ 5 (dashed-dotted line), T ¼ 50 (solid line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

g 0 ð0Þð1  cjÞ < 1; ð0; 0Þ is the only steady state. In the bistable regime, i.e. g 0 ð0Þð1  cjÞ > 1, three constant steady states exist, designated ð0; 0Þ; ðar ; br Þ; ðar ; br Þ. 2. (bs ¼ 0 for all Nontrivial Spatially-Dependent Solutions) For any 0 < Da , if as ðxÞ is spatially-dependent, bs ¼ 0 whenever c P 1 ~ for some 0 < k ~ < 1. and j P k, 3. (Existence of Nontrivial Spatially-Dependent Solutions) For any positive integer n, if Da < Dn2c , then for each s ¼ 1; 2; . . . n there exists a pair of steady-state solutions with s  1 interior extrema. 4. (Linear Stability of Homogeneous Steady States) (a) (Monostable regime) For g 0 ð0Þð1  cjÞ < 1; ð0; 0Þ is linearly stable for Da > Dc and linearly unstable for Da < Dc .

(b) (Bistable regime) For g 0 ð0Þð1  cjÞ > 1; ð0; 0Þ is linearly unstable for all 0 < Da . If, in addition, g 0 ðar  clbr Þ < 1, both ðar ; br Þ and ðar ; br Þ are linearly stable for all 0 < Da . If g 0 ðar  clbr Þ > 1, there exists a constant Dr , where 0 < Dr < Dc , such that both ðar ; br Þ; ðar ; br Þ are linearly stable for Da > Dr and linearly unstable for Da < Dr . 5. (Linear Stability of Spatially Dependent Steady States) For the solutions given in item 3, the monotone solutions (s ¼ 1) are the only solutions that can be linearly stable, and all others are linearly unstable.

5

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

a(x,t) 400

300

300

200

200

t

0

t

0.5 400

0.5 x

400

400

300

300

200

200

100 0 0

a,b(x,T)

0

0 0

1

t

t

0 0

0.5

100

100 −0.5

b(x,t)

0.5 x

1

−0.5 0

0.5 x

1

0.5

0

100 0.5 x

0 0

1

0.5 x

1

−0.5 0

0.5 x

1

Fig. 3. Each row represents the solution aðx; tÞ; bðx; tÞ using caa ¼ 10; l; c; j ¼ 1; q ¼ 0; Da ¼ 0:01; Db ¼ 1 in Eq. (1); the first two columns are the contour plots of aðx; tÞ; bðx; tÞ respectively, while in the last, aðx; TÞ (blue), bðx; TÞ (green) are shown for T ¼ 5 (dashed-dotted lines) and T ¼ 400 (solid lines); initial conditions in each row are random perturbations from the homogeneous equilibrium; while the transient behavior depends on the initial data, ultimately all solutions evolve to monotone patterns. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

b(x,t)

a(x,t)

200

150

150

100

100

t

0

t

0.5 200

50 0 0

0

0.5 x

0 0

1

0.5 x

1

−0.5 0

200

200

0.4

150

150

0.2

100

100

0

50

−0.2

50 0 0

a,b(x,T)

50

t

t

−0.5

0.5

0.5 x

1

0 0

0.5 x

1

−0.4 0

0.5 x

1

0.5 x

1

Fig. 4. Each row represents the solution aðx; tÞ; bðx; tÞÞ using caa ¼ 10; l; c; j ¼ 1; q ¼ 0; Da ¼ 0:01; Db ¼ 0:2 in Eq. (1); the first two columns are the contour plots of aðx; tÞ; bðx; tÞ respectively, while in the last, aðx; TÞ (blue), bðx; TÞ (green) are shown for T ¼ 5 (dashed-dotted lines) and T ¼ 200 (solid lines); initial conditions in each row are random perturbations from the homogeneous equilibrium; solutions evolve to both monotone and nonmonotone patterns, depending on initial data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The results of this theorem are illustrated in the following bifurcation diagrams shown in Fig. 5. Part 1 of the theorem was discussed in Section 1.1. 2.1. Preliminary results about the symmetry of steady-state solutions when b ¼ 0 It will be useful to introduce the change of variables u ¼ a  C, where C :¼ lcb. Then Eq. (4) becomes

0 ¼ Da u00  u þ gðuÞ  C; Z 1 cgðjuÞdx: C¼

ð5aÞ

Eðu; ux Þ :¼

Da 2 1 2 u  u þ GðuÞ  Cu ¼ c1 ; 2 x 2

where G is the antiderivative of g with Gð0Þ ¼ 0. We assume C and g 0 ð0Þ are such that there exist three constant equilibria to Eq. (5a), denoted eL ; ec ; eR with eL < ec < eR . When C ¼ 0; eL ¼ eR ; ec ¼ 0. We derive an elementary symmetry result in the following lemma. Lemma 2.2. For C ¼ 0, let uðxÞ be a solution of Eq. (5a), with boundary conditions u0 ð0Þ; u0 ð1Þ ¼ 0. Then

ð5bÞ

0

It is important to observe that Eq. (5a) is a Hamiltonian system with

ð6Þ

 u ¼ 0; and  either uðxÞ ¼ uð1  xÞ or uðxÞ ¼ uð1  xÞ.

6

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

0.5

a(0)

a(0)

0.5

0

0

0.5

0.5 0

20

40

60

80

0

20

40

1/Da

1/Da

(a)

(b)

60

80

a(0)

0.5

0

0.5 0

20

40 1/D

60

80

a

(c) Fig. 5. Bifurcation diagrams for steady-state solutions of Eq. (3) for parameter values g 0 ð0Þ ¼ 5=2, (a) c; j; l ¼ 1, (b) c; j ¼ 1; l ¼ 0:5, (c) c; j ¼ 1; l ¼ 0:10; (a) Represents the behavior in the monostable regime, where (0, 0) loses stability at Da ¼ 0:1520, and the monotone branch becomes stable; (b) Represents the first of two behaviors in the bistable regime, where the two nonzero steady states lose stability at Dr < Dc , which, from numerical results, corresponds to the point where the monotone branch becomes stable; (c) Represents the second possible behavior in the bistable regime, where the two nonzero steady states remain stable for all values of Da , but there exists a bifurcation along the monotone branch resulting in stable pattern formation.

TðaÞ ¼ Proof. In the case C ¼ 0, Eq. (6) simplifies to

Eðu; ux Þ :¼

Da 2 1 2 u  u þ GðuÞ ¼ c1 : 2 x 2

ð7Þ 0

We see that ðu; ux Þ ¼ ð0; 0Þ is a center, as long as g ð0Þ > 1. The solutions that satisfy the boundary conditions, u0 ð0Þ; u0 ð1Þ ¼ 0, rotate around the center either a 1/2 integer or whole integer number of times. All of these solutions satisfy u ¼ 0, as may be seen from the phase portrait, Fig. 6C. To prove the second part, we observe that the ðu; ux Þ phase plane associated to Eq. (7) has these symmetries. h We will use the following definitions

ð9Þ

Notice that Eq.(6) implies



dx duc

rffiffiffiffiffiffi Da 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 F c ðaÞ  F c ðuc Þ



¼

ð10Þ

Hence substituting Eq. (10) into Eq. (10)

TðaÞ ¼ 

Z

a

b

 dx duc : duc

ð11Þ

If, in addition uc ð1Þ ¼ b; uc ðxÞ is the monotone decreasing solution of the full boundary-value problem and

f ðuc Þ ¼ uc þ gðuc Þ; fc ðuc Þ ¼ f ðuc Þ  C; u2 Fðuc Þ ¼  c þ Gðuc Þ; 2

rffiffiffiffiffiffi Z a Da duc pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 b F c ðaÞ  F c ðuc Þ

ð8Þ

F c ðuc Þ ¼ Fðuc Þ  Cuc ; where G is the antiderivative of g with Gð0Þ ¼ 0. See Fig. 6A, B for graphs of fc ; F c . With these definitions, f is odd and has the same concavity properties as g. It is also useful to have the following time map, T, associated to Eq. (5a). Let uc ðxÞ be the solution to Eq. (5a) with uc ð0Þ ¼ a > ec ; u0c ð0Þ ¼ 0. There exists a unique value, b, on the same level set of the energy function as a, i.e. F c ðaÞ  F c ðbÞ ¼ 0. Then the map T is precisely the period of the orbit through a and b, which is defined as

TðaÞ ¼ 

Z

uc ð0Þ

uc ð1Þ



 Z 0 dx duc ¼  dx ¼ 1: du 1

ð12Þ

We can also use this integral constraint to find non-constant, nonmonotone solutions as well. For the solutions with a single interior extremum, the time map needs to take half as long, i.e. TðaÞ ¼ 12. Similarly,

TðaÞ ¼

1 n

will result in a solution with ðn  1Þ interior extrema. This information about the time map T will be used in the next section, in the proofs of parts 2 and 3 of Theorem 2.1.

7

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

B

fc,Fc

A 0.2 0.1 0 −0.1 −0.2

0.2 0.1 0 −0.1 −0.2

−0.5

uc’

C

0 u

0.5

−0.5

D

5

0

0

0.5

0

0.5

5

0

−5 −0.5

0 u

−5 −0.5

0.5

c

uc

E

0.5

F 0.5

0

0

−0.5

0

0.2

0.4

0.6

0.8

1

−0.5

0

0.2

0.4

0.6

0.8

1

x Fig. 6. A, B: Graphs of fc ; F c (solid blue, dashed black), see Eq. (8), for C ¼ 0; 0:05; C, D: Phase portraits for Eq. (5a) for C ¼ 0; 0:05; E,F:uc ðxÞ , the monotone decreasing solution to Eq. (5a), for C ¼ 0; 0:05. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2.2. Proof of parts 2 and 3 of Theorem 2.1 In this section, we prove parts 2 and 3 of Theorem 2.1, beginning with part 2. Let uc ðxÞ denote the monotone decreasing solution to Eq. (5a) with uc ð0Þ ¼ a and uc ð1Þ ¼ b. Throughout this section we keep C > 0 general. See Fig. 6 E,F for graphs of uc . We begin with a lemma, which will be used to by prove that the average value of uc ðxÞ is always positive. Lemma 2.3. uðxÞ þ uð1  xÞ > 0. Proof. a; b must lie on the same level set of the energy function, meaning F c ðaÞ  F c ðbÞ ¼ 0. Note that c

F c ðaÞ  F c ðe Þ > 0 and

F c ðaÞ  F c ðaÞ ¼ 2C a < 0:

sffiffiffiffiffiffi 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X ðxÞ ¼  F c ðaÞ  F c ðuc ðxÞÞ  F c ðaÞ  F c ðuc ð1  xÞÞ ; Da 0

1 X ðxÞ ¼ ðuc ðxÞ  gðuc ðxÞÞ þ uc ð1  xÞ  gðuc ð1  xÞÞ þ 2C Þ: Da 00

Consider

X00 ð0Þ ¼ u00c ð0Þ þ u00c ð1Þ ¼

1 ðða  bÞ  gðaÞ þ gðbÞ þ 2C Þ Da

and let

Hða; bÞ ¼ ða  bÞ  gðaÞ þ gðbÞ þ 2C ¼ f ðaÞ þ f ðbÞ þ 2C: From the constraint that F c ðaÞ  F c ðbÞ ¼ 0, we get

FðaÞ  FðbÞ ¼ C: aþb Thus,

As F c ðuÞ is decreasing for eL < u < ec ; b is the unique zero of F c ðaÞ  F c ðuÞ for this range. Therefore, b must lie between a and ec , or a  b > 0. Define

Hða; bÞ ¼ f ðaÞ þ f ðbÞ þ 2

XðxÞ :¼ uc ðxÞ þ uc ð1  xÞ:

FðaÞ  FðbÞ ¼

XðxÞ is clearly even about x ¼ 1=2, by definition, with Xð0Þ ¼ a  b > 0. We will show X has local minima at the endpoints and is positive for all x 2 ½0; 1. We assume b < 0. Otherwise, uc ðxÞ > 0 for all x 2 ½0; 1 and X > 0 by assumption. Direct computation gives

FðaÞ  FðbÞ : aþb

As f is concave down for positive values,

Z b

a

f ðuÞdu P

Z b

a

f ðaÞ  f ðbÞ ðu  aÞ þ f ðaÞdu ab

f ðaÞ  f ðbÞ ðb  aÞ2 þ f ðaÞða  bÞ ab 2 ða  bÞ ¼ ðf ðaÞ þ f ðbÞÞ: 2 ¼

ð13Þ

8

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

Therefore,

Using the results of the average value of steady-state solutions, we now address the question of the number of such solutions. We will use the relationship between a and Da given by Eq. (9), (12) to answer these questions. As C ¼ 0; b ¼ a. Using the change of variables, u ¼ a sin h, we can rewrite TðaÞ as

FðaÞ  FðbÞ aþb ab P f ðaÞ þ f ðbÞ þ ðf ðaÞ þ f ðbÞÞ aþb 2 ðbf ðaÞ þ af ðbÞÞ; ¼ aþb

Hða; bÞ ¼ f ðaÞ þ f ðbÞ þ 2

rffiffiffiffiffiffi Z Da p=2 a cos h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dh 2 p=2 FðaÞ  Fða sin hÞ pffiffiffiffiffiffiffiffiffi Z p=2 a cos h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dh: ¼ 2Da FðaÞ  Fða sin hÞ 0

TðaÞ ¼

which is positive for a > b > 0, giving X00 ð0Þ > 0. Suppose Xðx Þ ¼ 0 for some x 2 ð0; 1=2Þ and is the smallest such point. Since Xð0Þ > 0; X0 ð0Þ ¼ 0; X00 ð0Þ > 0, then there must exist ^ x 2 ð0; x Þ so that Xð^ xÞ > 0; X0 ð^ xÞ ¼ 0, and X00 ð^ xÞ < 0. Let ^ ^ v :¼ uc ðxÞ and w :¼ uc ð1  xÞ. Then, X0 ð^xÞ ¼ 0 implies F c ðv Þ  F c ðwÞ ¼ 0. Yet, X00 ð^ xÞ ¼ f ðv Þ þ f ðwÞ þ 2C > 0 from above argument, a contradiction. Hence, XðxÞ > 0. h

As in [61], we can use this expression to obtain T 0 ðaÞ

T 0 ðaÞ ¼

Using Lemma 2.3 above, we prove the average value of uc ðxÞ is always positive. Lemma 2.4.

uc :¼

Z

uc ðxÞdx > 0:

uc ¼

uc ðxÞdx þ

Z

0

¼

Z

1

uc ðxÞdx ¼

Z

1=2

0

1=2

uc ðxÞdx þ

Z

0

uc ð1  xÞdx

1=2

ðFðaÞ  FðuÞÞ3=2

0

du:

pffiffiffiffiffiffi p Da ffi: limTðaÞ ¼ pffiffiffiffiffiffiffiffiffiffiffi a!0 F 00 ð0Þ 00

1=2

uc ðxÞ þ uc ð1  xÞdx:

0

From Lemma 2.3, the integrand is always positive, hence uc > 0. h Next, we prove a result regarding the average value of gðuc ðxÞÞ  cgðjuc ðxÞÞ, as this is equivalent to uc for solutions of the full system Eq. (5). ~ 2 ð0; 1Þ such that for c P 1 and j P k, ~ Theorem 2.5. There exists a k

Z

a

Corollary 2.6. T is an increasing function in a. In addition, the limit as a ! 0 can be calculated, yielding

Proof. Consider the average of uc , 1=2

 cos hdh pffiffiffiffiffiffiffiffiffi Z a 2Da FðaÞ  FðuÞ  a=2f ðaÞ þ u=2f ðuÞ

The numerator of the integrand is always positive, again as it is an immediate consequence of the concavity of f. Hence we conclude T is a strictly increasing function in a.

1

0

Z

¼

pffiffiffiffiffiffiffiffiffi Z p=2 FðaÞ  Fða sin hÞ  a=2ðF 0 ðaÞ  F 0 ða sin hÞ sin hÞ 2Da ðFðaÞ  Fða sin hÞÞ3=2 0

0

1þg ð0Þ As the monotone solution satisfies TðaÞ ¼ 1; Dc :¼ F pð0Þ cor2 ¼ p2 responds to a pitchfork bifurcation point, where the curve of monotone solutions emerges. Therefore, for any 0 < Da < Dc , the monotonicity of T gives a unique monotone increasing and decreasing solution. In addition, we can consider the existence of the solution which is monotone on x 2 ½0; 1=n, i.e. must solve TðaÞ ¼ 1=n. Therefore, for 0 < Da < Dc =n2 , this solution will exist. For any fixed Da , a finite number of paired solutions exists, with the number depending on the smallest n such that Da < Dc =n2 . This completes the proofs of parts 2 and 3 of Theorem 2.1.

1

½gðuc ðxÞÞ  cgðjuc ðxÞÞdx < 0:

0

Proof. Let

KðkÞ ¼

Z

1

gðuc ðxÞÞ  cgðjuc ðxÞÞdx:

0

R1 R1 Kð0Þ ¼ 0 gðuc ðxÞÞ > C > 0 and Kð1Þ ¼ ð1  cÞ 0 gðuc ðxÞÞ < 0 from ~ 2 ð0; 1Þ such that KðkÞ ~ ¼ 0. Lemma 2.4. Therefore there exists a k In addition, for any j > 1, Z 1 Z 1 gðjuc ðxÞÞdx > gðuc ðxÞÞdx; ð14Þ 0

0

~ and c P 1. h so that KðjÞ < 0 for all j > k The proof of part 2 of Theorem 2.1 follows as an immediate corollary. Proof. For a solution of Eq. (5), we have uc > 0, but for c P 1 and ~ k P k,

uc ¼

Z

gðuc ðxÞÞ  cgðjuc ðxÞÞdx < 0;

a contradiction. Hence, there is no solution satisfying Eq. (5) for C > 0. Similarly, no solution can exist for C < 0. h

2.3. Proof of parts 4 and 5 of Theorem 2.1 We will begin by analyzing the stability of the constant steadystate solution. In the monostable regime, (0, 0) is the only constant steady state. Linearizing the system around (0, 0), we have eigenvalue/eigenfunction k; ðw; v Þ satisfying

kw ¼ Da wxx  w þ g 0 ð0Þðw  lcv Þ; Z 1 kv ¼ lv þ jg 0 ð0Þ ðw  lcv Þdx:

ð15Þ

0

v is spatially independent, jg 0 ð0Þ v¼w : k þ l þ jlcg 0 ð0Þ As

Therefore, ðcosðpnxÞ; 0Þ is an eigenfunction corresponding to eigenvalue kn ¼ Da p2 n2  1 þ g 0 ð0Þ for n ¼ 1; 2; 3; . . . There also exist two eigenvalues, k0;1 ; k0;2 associated to the spatially constant eigenvectors, which both have negative real parts by the monostability assumption. Taking Da as a bifurcation parameter, consider decreasing values of Da :  k0;1 ; k0;2 < 0 for all Da . 0 ð0Þ  Dc :¼ 1þg is the first simple bifurcation point, in that p2 k1 ¼ 0; ki < 0; i ¼ 2; 3; . . . when Da ¼ Dc , and k1 > 0 for Da < Dc .

9

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

 When Da ¼ D22c ; k2 ¼ 0; ki < 0; i ¼ 3; 4; . . ., and k2 > 0 for Da < Dc =22 .  Similarly, Dn2c is the n  th simple bifurcation point, where kn ¼ 0 and will be positive for all smaller Da .

In this section, the shape of the monotone decreasing solutions of Eq. (4) is briefly discussed. Two classic pattern types are seen in such reaction–diffusion systems, often referred to as ‘‘spikes’’ and ‘‘plateaus’’, with definitions given in [62]. The sign of the fourth derivative determines the type, with positive values associated to spikes, negative to plateaus. Fig. 7(a) gives the value of að4Þ ð0Þ for the monotone decreasing solutions of Eq. (4) as a function of Da , with aðxÞ=að0Þ shown in (b). að4Þ ð0Þ < 0 for 0 < Da < 0:1 and attains a minimum near Da ¼ 0:02. Thus, for Da ¼ 0:001 and 0.02, the associated front solutions are both plateaus. For larger values of Da , the front solutions become spikes, as they resemble the corresponding eigenmode, cosðpxÞ, near the bifurcation point, Dc ¼ 0:1520. However the stability of the front solutions does not change with the type.

3. Linear analysis of full coupled system We now consider the full coupled system, Eq. (1), again in the absence of f. While a complete picture of all steady states and their

10

1

0

0.8

10

0.6

20

0.4

30

0.2

D = 0.001 a

a/a(0)

(4)

a (0)

Thus, for Da > Dc ; ð0; 0Þ is stable and is unstable for Da < Dc . This summarizes the stability of the only constant solution in the monostable regime, completing the proof of part 4(a). However, by virtue of the reduced system being a shadow system, only constant and monotone solutions can be stable steady states [59], finishing part 5. Numeric results support the stability of the monotone branch. In the bistable regime, (0, 0) is a steady state, but always unstable by assumption. For the nonzero constant solutions, ðar ; br Þ; ðar ; br Þ; k0;1 ; k0;2 < 0 by assumption. If 1 þ g 0 ðar  lcbr Þ > 0, then there exists a Dr < Dc so that k1 > 0 when Da < Dr , completing part 4(b). Numerical results indicate when such a bifurcation occurs, it coincides with a bifurcation along the monotone branch, stabilizing such solutions. However, even when the constant solutions remain stable for all Da , bifurcation still occurs along the monotone branch, leading to three stable steady-state solutions, see Fig. 5.

2.4. Shape of solutions

40

=0.02 =0.12

0 0.2

50

0.4

60

0.6

70

0.8

80 0

0.02

0.04

0.06

0.08

Da

(a)

0.1

0.12

0.14

1 0

0.2

0.4

0.6

0.8

1

x

(b)

Fig. 7. (a) að4Þ ð0Þ for the monotone decreasing solution of Eq. (4) with l; c; j ¼ 1; g 0 ð0Þ ¼ 2:5; (b) a=að0Þ, the scaled front solutions for Da ¼ 0:001; 0:02; 0:12; For Da ¼ 0:001; 0:02; að4Þ ð0Þ < 0, indicating a plateau solution, while for Da ¼ 0:12; að4Þ ð0Þ >, indicating a spike.

(a)

(b)

Fig. 8. (a) Characteristic wavelength for linear system associated to Eq. (1) with c; j; l ¼ 1; g 0 ð0Þ ¼ 10; q ¼ 0; Black lines designate bifurcation curves, kn ¼ 0, see Eq. (17), and color designates the maximum eigenvalue; (b) Regions of stability for stable monotone (SM) and stable non-monotone (SNM) solutions.

10

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

stabilities cannot be computed as in the shadow system, linear stability analysis of the constant solution provides evidence for the long term dynamics. We will assume Db > Da . To each eigenfunction cosðpnxÞ, there correspond two eigenvalues, kn;1 ; kn;2 , given explicitly by

  pffiffiffi  p2 n2 ðDa þ Db Þ þ 1 þ l þ g 0 ð0Þðljc  1Þ þ d ; kn;1 ¼ 2   pffiffiffi  p2 n2 ðDa þ Db Þ þ 1 þ l þ g 0 ð0Þðljc  1Þ  d ; kn;2 ¼ 2

2 2 d ¼ ðDb  Da Þp2 n2 þ l  1 þ ðcjl þ 1Þg 0 ð0Þ  4cjl½g 0 ð0Þ : ð16Þ As kn;1 will have the largest real part, we will let

kn :¼ kn;1 :

Parameter set

x0 Mean

Stdev

Max

Min

Range

1 2 3 4 5 6 7 8 9 10 11

0.5561 0.5611 0.5215 0.5561 0.5561 0.5314 0.5891 0.5858 0.5314 0.5561 0.5512

0.0565 0.0566 0.0682 0.0565 0.0565 0.0311 0.0856 0.0895 0.0311 0.0565 0.0552

0.6337 0.6436 0.6040 0.6337 0.6337 0.5743 0.7030 0.7030 0.5743 0.6337 0.6238

0.4950 0.4950 0.4257 0.4950 0.4950 0.4950 0.4950 0.4851 0.4950 0.4950 0.4951

0.1387 0.1486 0.1783 0.1387 0.1387 0.0793 0.2080 0.2179 0.0793 0.1387 0.1287

ð17Þ

Straightforward computation reveals all eigenvalues are real for a range of parameters. Lemma 3.1. kn;1 ; kn;2 is real for 1. 2.

Table 1 aðx; 200Þ was computed numerically for Eq. (18) for eleven sets of parameter values; x0 designates the inflection point for aðx; 200Þ, giving a definition for the location of the front; for each parameter set, mean, standard deviation, max, min, and range of x0 was computed for six functions f(x) depicted in Fig. 10(a).

l P 1, or l l < 1 and g 0 ð0Þ > ð11 pffiffiffiffiffiffi 2 ; cjl – 1, cjlÞ

For fixed values of c; j; l; g 0 ð0Þ, the bifurcation curves, i.e. kn ¼ 0, are easily computed and are plotted in the Da ; Db plane as black curves in Fig. 8(a). Any point lying above the curve corresponds to diffusion values which produce a positive nth eigenvalue for the linear system. This gives the number and type of unstable eigenmodes. The eigenmode associated with the largest positive eigenvalue is also important, since it yields the characteristic wavelength. The largest eigenvalues were computed numerically and designated by colors in Fig. 8(a). The dark blue region represents values of Da ; Db with k0 the largest, medium blue region k1 the largest, and so on. For a fixed value of Db , as Da decreases, more eigenvalues become positive and the characteristic wavelength increases. However, for a fixed value of Da , increasing Db decreases the characteristic wavelength. Thus, for large enough Db ; k1 is the largest eigenvalue, corresponding to the monotone eigenmode. Numeric simulations support the linear stability results in predicting the number and stability of the nonconstant steady states. Fig. 8(b) demonstrates the regions in which stable monotone (SM) and stable non-monotone (SNM) solutions exist. Solutions were considered stable if patterns persisted past T ¼ 800 with noise introduced into the numeric solver. As in, Figs. 3, 4, decreasing Db can stabilize nonmonotone steady-state patterns. But for large values of Db , the characteristic wavelength is one and only the monotone solutions are stable, the same as that for the shadow system. 4. Robustness We now consider the behavior of the systems, Eqs. (1), (3), with variable, spatially dependent functions, f(x). The introduction of functional inputs, f(x), into the sigmoidal production term accounts for known biological regulators that vary over space, but are relatively fixed over time. As Bcd gradients have been approximated by exponential functions [37,63], we use noisy exponentials for f(x). In both the shadow and full system, the monotone steady-state solution remains stable and near the solution in the absence of f. Of particular importance in creating robust patterning under variability of f(x) is the magnitude of c, which measures the strength of the inhibitory interactions in the network. Similar robustness was

1 2 3 4 5 6 7 8 9 10 11

g0 (0)

l

c

q

j

2.5 5 1.25 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5

1 1 1 2 0.5 1 1 1 1 1 1

1 1 1 1 1 2 0.5 1 1 1 1

1 1 1 1 1 1 1 2 0.5 1 1

1 1 1 1 1 1 1 1 1 2 0.5

demonstrated in [26], but without discussion on the primary parameters responsible for the property.

4.1. Shadow system robustness Without additional input functions f to the synthesis production g, the symmetry of g forces symmetry of the steady states around the midpoint, x ¼ 1=2, see Section 2.1. We define x0 as the inflection point of the monotone decreasing, steady-state front. With the addition of nontrivial f to Eq. (3), we consider

at ¼ Da axx  a þ gða  lcb þ qf ðxÞÞ; Z 1 b0 ¼ lb þ gðjða  lcb þ qf ðxÞÞdx:

ð18Þ

0

x0 was computed numerically for a range of parameters using f ðxÞ ¼ Mekx þ gðxÞ, for ðM; kÞ ¼ ð0; 0Þ; ð1; 3Þ; ð1; 4Þ; ð1; 5Þ; ð1; 10Þ; ð0:5; 5Þ; gðxÞ a random number between ð0:1; 0:1Þ for each x. In doing so, the effect of individual parameters on the robustness of the pattern can be identified. Results are listed in Table 1. Parameter sets 2–11 are viewed as deviations from the reference set 1. x0 remained robust for all parameter sets, i.e. the monotone steady-state front exists and remains stable under changes to f(x). The only parameter changes that reduced the standard deviation and range of x0 were increasing c (set 6) and decreasing q (set 9). Clearly, the smaller q is compared to the other terms, the less effect f(x) has relative to other interactions. When q ¼ 0; f ðxÞ has no effect on the dynamics, thus this is viewed as a trivial case. However, the larger c, the larger the inhibitory interactions in the network. We speculate that with increasing negative feedback, the network can filter larger variations of f(x). Further increasing c to 10, the mean and standard deviation of x0 are 0.4918, 0.0253, respectively. The changes to l; j had little effect on x0 , when it exists. Fig. 10 demonstrates the robust patterning seen when c ¼ 5. Fig. 9 shows the decrease in standard deviation of x0 when increasing values of c for four parameter sets.

11

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

0.09

/c ðxÞ ¼ wðxÞ  cðx  1=2Þ;

0.08

pffiffiffiffiffiffi where ca ¼ 1= Da . With continuity and differentiability of a at x0 ,

ð22Þ

aðx0 Þ ¼ wðx0 Þ;

0.07

thus Eq. (20) becomes

0.06

Stdev

/c ðx0 Þ þ f ðx0 Þ ¼ 0:

0.02

For details, see Appendix B. The graphs of w; /c are depicted in Fig. 11 for various values of Da . For very small values of Da ; w  0 away from x ¼ 0 and x ¼ 1. Similarly, /c  cðx  1=2Þ away from the endpoints. For a linear input function, f ðxÞ ¼ 1=2  mx for some constant m; 1=2  x0 is approximated as

0.01

1=2  x0 

0.04 0.03

0

1

2

3

4

5

γ

6

7

8

9

10

Fig. 9. Standard deviation of x0 as a function of c for Eq. (18); Da ¼ 0:01, f depicted in Fig. 10(a); l; j; q ¼ 1; g 0 ð0Þ ¼ 2:5 (blue, solid, circles), l; j ¼ 1; q ¼ 2; g 0 ð0Þ ¼ 2:5 (green dashed +), l ¼ 2; j ¼ 0:5; q ¼ 1; g0 ð0Þ ¼ 5 (red dotted star), l ¼ 0:5; j ¼ 2; q ¼ 2; g0 ð0Þ ¼ 2:5 (cyan dashed-dotted diamond). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Now we obtain an explicit characterization of the front location when taking g ¼ H  1=2, with H the Heaviside function. The equations governing the front position, x0 , are

aðx0 Þ  lcb þ qf ðx0 Þ ¼ 0 Z 1 Hðjða  lcb þ qf ðxÞÞdx  1=2 ¼ 0:  lb þ

ð19Þ

0

As we are assuming a monotone decreasing front solution for a,

b¼

ð23Þ

0.05

1

l

ðx0  1=2Þ;

thus Eq. (19) becomes

aðx0 Þ  cðx0  1=2Þ þ qf ðx0 Þ ¼ 0:

ð20Þ

Define

wðxÞ ¼

sinhðxca Þ coshðð1  xÞca Þ 1  sinhðca Þ 2

ð21Þ

and

m1 : 2ðc þ mÞ

Similarly, for an exponential input function, f ðxÞ ¼ ekx for some constant k, we have the approximation

1=2  x0  

1

cek=2 þ k

:

As j1=2  x0 j represents the shift of the front location away from the midpoint, for large c; j1=2  x0 j  1. Therefore, with stronger inhibition, b can filter out the positional variability induced by f(x) on the final pattern for a. While this was assumed for g ¼ H  1=2, numeric results indicate close correlation for smooth g as well. 4.2. Robustness to parameter changes In this section, we explore robustness in the general model. Recall Eq. (1) originated from the more general model, Eq. (A.2) given in Appendix A. By selecting random parameter values in Eq. (A.2), we first demonstrate that robustness occurs in less than 1% of such equations. In contrast, with the addition of conditions Eqs. (B.3) and (B.6) obtained in Appendix B, robustness is seen in nearly one-third of all such systems. For the first robustness test, Da ¼ 0:01 and g 0 ð0Þ ¼ 2:5 are fixed, while cij ; l; ri ; j are all chosen randomly in the range of ð5; 5Þ, except cfa ; caa ; l 2 ð0; 5Þ. When f ¼ 0, the front solution exists in 9 of the 1000 systems. Of these 9, in only 5 did x0 remain within (0.4, 0.6) when f ¼ e5x . Thus, only 0.5% of sets demonstrated the desired robustness, Fig. 12(a). On the other hand, when cfa ; caa ; l; j; ri are again chosen randomly in same ranges, but the remaining parameters are determined by Eqs. (B.3) and (B.6), the front solution exists and is stable in over 300 sets. In addition, 0.5

1

0.4 0.3

0.8

0.2 0.1

f

a

0.6 0.4

0 0.1 0.2

0.2

0.3 0.4

0 0

0.2

0.4

0.6

0.8

1

0.5

0

0.2

0.4

0.6

x

x

(a)

(b)

0.8

1

Fig. 10. Numeric simulation of steady-state solutions to Eq.(18), Da ¼ 0:01; c ¼ 5; j; l; q ¼ 1; g 0 ð0Þ ¼ 2:5; f ðxÞ and aðx; 200Þ; Even with large variability in f, the monotone solution remains the only stable steady state, and the location of the front is robust.

12

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

0.5

0.5

D = a

1

0.01 0.05 0.1 0.5 1

0

0.5

0

0.2

0.4

x0

0.6

0.8

0

0.5

1

0

0.2

0.4

(a)

0.6

x0

0.8

1

(b)

Fig. 11. Graph of (a) wðx0 Þ, (b) /1 ðx0 Þ for various choices of Da ; For small values of Da ; wðx0 Þ is well approximated by zero away from the endpoints, while /1 ðx0 Þ is well approximated by ð1=2  x0 Þ.

4

f=0 f=e

450 400

5x

350

Frequency

3

Frequency

f=0 f=e

5x

2

300 250 200 150

1

100 50

0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

0

1

0

0.1

0.2

0.3

0.4

0

0.5 x

0.6

0.7

0.8

0.9

1

0

(a)

(b)

Fig. 12. Stacked bar graphs for the frequency, per 1000 parameter sets, of the position of the front, x0 , for f ¼ 0 (blue) and, of those sets with x0 2 ð0; 1Þ for f ¼ 0; x0 for f ¼ e5x (grey) in Eq.(18); Da ¼ 0:01; g 0 ð0Þ ¼ 2:5; (a) cij ; l; r i ; j are all chosen randomly, cfa ; caa ; l 2 ð0; 5Þ, all others 2 ð5; 5Þ; (b) cfa ; caa ; l; j; ri chosen randomly in ranges given in (a), but remaining parameters are determined by Eqs. (B.3) (B.6). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

0.5 1

0.4 0.3

0.8

0.2 0.1

f

a

0.6 0.4

0 0.1 0.2

0.2

0.3 0.4

0 0

0.2

0.4

0.6

0.8

1

0.5

0

0.1

0.2

0.3

0.4

0.5

x

x

(a)

(b)

0.6

0.7

0.8

0.9

1

Fig. 13. Numeric simulation of steady-state solutions to Eq.(1), Da ¼ 0:01; Db ¼ 1; c ¼ 5; j; l; q ¼ 1; g 0 ð0Þ ¼ 2:5; f ðxÞ and aðx; 200Þ; When the monotone solution is stable, the position of the front is largely robust to variability in f, but unlike Fig. 10, other nonmonotone patterns can exist and be stable; (a) same as Fig. 10(a).

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

or equivalently,

0.12

wðx0 Þ þ qf ðx0 Þ þ ra ¼ 0; 0.1

ð24Þ

0.06

where w is defined in Eq. (21). Suppose x0 ¼ 1=2 is a solution to Eq. (24). As wð1=2Þ ¼ 0; w0 ð1=2Þ ¼ ca = sinhðca Þ; w becomes very flat around x0 ¼ 1=2 for large ca (i.e. small Da ), see Fig. 11. Therefore, even small deviations of f can result in large deviation of x0 from 1/2 if f ð1=2Þ does not remain fixed. Thus the dynamics of the scalar equation are not robust to changes in f.

0.04

5. Conclusions

0.08 Stdev

13

0.02

0

0.2

1

Db

10

100

Fig. 14. Standard deviation of x0 as a function of Db for Eq.(1); l; j; q ¼ 1; g 0 ð0Þ ¼ 2:5, f depicted in Fig. 10(a); c ¼ 1; Da ¼ 0:01 (blue, solid, circles); c ¼ 1; Da ¼ 0:025 (green dashed x); c ¼ 1; Da ¼ 0:05 (red dotted star); c ¼ 5; Da ¼ 0:01 (cyan dashed-dotted triangle), no point exists at Db ¼ 0:2 due to the instability/non-existence of front solution at that value. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

261 demonstrated robustness to the change in f. Thus, the relationships between parameters given by Eqs. (B.3) and (B.6) largely govern the location and robustness of solutions, not the exact values of the parameters. 4.3. Robustness in full system For large values of Db , numerics and linear stability results indicate the dynamics for the full system resemble that of the shadow system. This remains true in the presence of functions, f(x). For Da ¼ 0:01; Db ¼ 1, values of x0 were numerically computed, with near identical results to that in Table 1. However, as Db decreases, multiple stable steady-state patterns may exist, indicated by Fig. 13, which is also seen with nontrivial f. Fig. 13 again depicts the steady-state profiles under noisy f ; c ¼ 5. As opposed to Fig. 10, solutions can converge to nonmonotone patterns, but still show reduced variability in front location with increasing values of c. Fig. 14 demonstrates two trends regarding the standard deviation with respect to the diffusion rates, Da ; Db . First, the standard deviation of x0 decreases with increasing values of Db . Thus, the shadow system, which is obtained as the limit when Db ! 1, is the most robust. Second, the first three parameter sets show a decrease in standard deviation with increasing values of Da . While a small set, these parameter sets are representative of the behavior of the system in general. 4.4. Second protein necessary for robustness We will now show that the dynamics for a scalar equation cannot be robust to changes in f(x). As in Section 4.1, we will consider g ¼ H  1=2, to obtain explicit formula for x0 , the front border. Consider the steady-state scalar equation for the dynamics of a single protein, aðxÞ,

0 ¼ Da axx  a þ gða þ qf ðxÞ þ r a Þ: Thus x0 must satisfy

aðx0 Þ þ qf ðx0 Þ þ r a ¼ 0

We have provided a model exhibiting three key characteristics of the Bcd–Hb system, in that solutions converge to a monotone pattern symmetric around the midpoint of space, despite changes in (a) spatial scalings, (b) initial conditions, and (c) input f, within a class of functions. Bcd profiles are represented by f(x), given as exponentials fixed in time, a, representing the anterior Hb domain, and b, a second protein. By constructing a ‘‘minimal model’’ capable of such robust pattern formation, we find  A two component network is the smallest such network exhibiting the robust properties, with topology indicated by Fig. 1.  The magnitude of the inhibitory interactions primarily controls the buffering capacity of the system against perturbations of f(x).  Increasing the diffusion rate for the second protein increases the robustness of the system and destabilizes non-monotone solutions. Under these conditions, a is robust to changes in positional information and produces a precise pattern, not dictated by a simple threshold effect. Fig. 13 demonstrates the capacity of such a model in filtering noisy input. Each of the above points is addressed in more detail, as follows. 5.1. b and Network topology As indicated by Section 4.4, the second protein, b, is critical for the robust features. In this model, b can be spatially homogeneous (as in the shadow system, Eq. (3)), form an anterior-to-posterior gradient (as in Eq. (1), Network Type I), or form a posterior-toanterior gradient (as in Eq. (1), Network Type III, data not shown). The consistent interaction coefficients in both Network Types I and III indicate self-inbition of b and an opposite interaction with a, i.e. cAB cBA < 0 and cBB < 0, see Fig. 1. With regard to Bcd–Hb system, current literature still cites the reduced positional variability of Hb and other downstream genes, although with debate. One such model proposed requires a posterior gradient, forming a mirror image to that of Bcd [35–37]. The two oppositely directed gradients are balanced at the midpoint of the embryo, indicating where Hb should drop. While such a system could account for the Hb border at the midpoint of the embryo, it heavily relies on the symmetry of the two gradients, as well as correlated perturbations [64,26]. In addition, known candidates for the posterior gradient have been eliminated, and no other such gradient has been identified [15]. Similarly, models requiring active transport [64] have not been experimentally verified. If b represents a single factor in the biological system, the most likely candidates are the zygotic genes Kr (Krüppel) and kni (knirps), both known repressors of hb expression [65–67]. There is evidence to suggest auto-repression of Kr and kni [68,69], as well as Kr activation by Hb at low levels [70,13,29], which would be consistent with the network topology of the model. However, indi-

14

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

vidual mutations to Kr and kni did not significantly alter location or standard deviation of the Hb border [15], but it was disrupted in double mutants [22]. More likely, the role of b may represent a subnetwork in the biological system, such as Kr-kni, cumulatively acting as an inhibitor of hb. Alternatively, as b can be spatially homogeneous, it could be overlooked as a contributing factor for the Hb dynamics. 5.2. Inhibitory strength, c From Section 4, results indicate that increasing the parameter c increased the robustness of the location of the front to changes in the function, f(x) (Table 1). The magnitude of c is proportional to the magnitude of the inhibitory coefficients in the network, as cba ¼ caa lc and cbb ¼ jcba . In terms of dimensional parameters,

c¼

 cBA RA kB : cAA kA RB

There are two natural interpretations for increasing the value of c. BA j First, c increases as jccAA increases, i.e. as the inhibitory affect of B on A increases relative to the auto-activation rate of A. Second, as Ri =ki is the saturation level for protein i; c also increases as the ratio of saturation levels for A to B increases. Thus, the model predicts increasing robustness with increasing inhibition or increasing the amount of A compared to B. The importance of the inhibitory interactions is supported by biological evidence, which indicates strong, mutual repression of gap genes [66,71]. In addition, the interpretation of inhibitory strength as the driving force in reduced Hb variability is in alignment with the results of [22], where the balancing of inhibition (Krüppel and Knirps) and activation (Bcd) determines to location for the Hb border. 5.3. Diffusion and stability Recall that the shadow system, Eq. (3), was obtained from Eq.(1) by considering instantaneous diffusion of b. The non-dimensional diffusion rate, Da ¼ DA =ðkA L2 Þ, and thus changes in Da are equivalent to spatial scalings of the original dimensional variable. Thus, with regard to robustness to spatial scaling, as long as Da < Dc , the first bifurcation point, the symmetric monotone solutions exist and are the only stable steady states (Fig. 5). The only dependence on initial conditions is whether solutions will converge to the monotone decreasing or increasing solution. As a consequence of the symmetry of solutions, initial conditions are equally divided into the basin of attraction for each, resulting in a large class of functions with the same long-term dynamics. Non-constant, nonmonotone steady states exist for Da < Dc =n2 , but are never stable, and thus may contribute to transient behavior of the system only. Recent studies have also indicated a reduction of variability of gene patterns over time, termed ‘‘canalization’’, leading support to considering developmental patterns as steady states of dynamical systems [72,21]. In the full system, the location and robustness for the monotone front is preserved, but multiple stable solutions can exist with more comparable rates of diffusion (Fig. 13). Numerics predict the stability of nonmonotone steady-state patterns, in contrast to the dynamics of the shadow system. However, linear stability analysis still predicts with large values of Db , the characteristic wavelength is one, corresponding to the monotone eigenmode (Fig. 8). Therefore, with Db sufficiently large but finite, the dynamics resemble that of the shadow system. While results in Section 2 focused on the shadow system dynamics, Fig. 14 indicates that increasing Db increases the robust-

ness of the system. Thus, the shadow system is the ‘‘most robust’’, but similar properties are found for systems with finite values of Db . 5.4. Final remarks While undoubtedly a simplification, this model for the Bcd–Hb dynamic has identified a two-component system able to replicate all desired robust features. A second, dynamic protein is essential, where the strength of its inhibition on Hb allows the system to effectively filter noisy positional input, such as Bcd and other known regulators. With faster rates of diffusion for the second protein, robustness is increased, supporting the reduction to the shadow system, which also acts to destabilize nonmonotone solutions. With these features, the long term dynamics act almost independently of initial conditions and inputs. Such a simple system could serve as a building block to better understand the complexity of the true mechanism for positional precision under noisy conditions. Acknowledgments I would like to thank R. Forman, my thesis advisor who supported and directed the initial project. I would also like to thank the helpful discussions and advice from T. Kaper, B. Mukhopadhyay, J. Reinitz, and reviewers. This work was partially supported by the Center for BioDynamics at Boston University and the NSF DMS 0602204 (EMSW21-RTG). Appendix A. General model and relation Let AðX; TÞ; BðX; TÞ represent the concentrations of the two proteins, A and B, at time T and location X. We assume zero flux condition AX ; BX ¼ 0 when X ¼ 0; L. We are modeling the developmental time window between cell cycles 11 and 14, before cellularization, so proteins diffuse freely, where DA ; DB are constants determining the rate of diffusion of each protein. We assume a first order decay of each protein, where kA ; kB are the rate constants. The production of A is a function of the form RA gðuA Þ, where RA is the maximal production rate of A; gðuA Þ is monotonically nondecreasing, sigmoid function and uA is a linear combination of A; B. Specifically, uA ¼ cAA A þ cBA B þ r a þ cFA FðXÞ, where cAA ; cBA denote the pairwise influences of protein A on the production of protein A and the influence of protein B on the production of protein A, respectively. Uniform maternal input is represented by ra , while maternal time-independent gradients are represented by FðXÞ, with cFA denoting the pairwise effect of F on the production of A. Therefore, the general model is of the form

@A @2A ¼ DA 2  kA A þ RA gðcAA A þ cBA B þ cFA FðXÞ þ ra Þ; @T @X @B @2B ¼ DB 2  kB B þ RB gðC AB A þ cBB B þ cFB FðXÞ þ r b Þ; @T @X 0 6 X 6 L;

ðA:1Þ

T P 0;

AjT¼0 ¼ A0 ðXÞ; BjT¼0 ¼ B0 ðXÞ;

@A

@B

¼ ¼ 0: @X X¼0;L @X X¼0;L With the following change of variables,

A ¼ aRA =kA ;

T ¼ t=kA ;

we obtain the system

X ¼ xL;

B ¼ bRB =kA ;

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

@a @2a ¼ Da 2  a þ gðcaa a þ cba b þ cfa f ðxÞ þ r a Þ; @t @x @b @2b ¼ Db 2  lb þ gðcab a þ cbb b þ cfb f ðxÞ þ r b Þ; @t @x 0 6 x 6 1; t P 0 ajt¼0 ¼ a0 ðxÞ; bjt¼0 ¼ b0 ðxÞ;

@a

@b

¼ ¼ 0; @x @x x¼0;1

Uðx0 ; y0 Þ þ r a =caa þ cfa =caa f ðx0 Þ ¼ 0; Wðx0 ; y0 Þ þ r b =cab þ cfb =cab f ðy0 Þ ¼ 0: ðA:2Þ

cFA ¼ cfa ;

cFB ¼ cfb ;

caa ¼ cAA RA =kA ;

l ¼ kB =kA ;

f ðxÞ ¼ FðX=LÞ;

cab ¼ cAB RA =kA ;

cba ¼ cBA RB =kA ;

Using first order linear approximations for U; W around ð1=2; 1=2Þ, with

@U c coshðcb =2Þ2 ð1=2; 1=2Þ ¼ cba =caa b ; @y0 sinhðcb Þ @W ca ð1=2; 1=2Þ ¼ coshðca =2Þ2 ; @x0 sinhðca Þ

where

Db ¼ DB =ðkA L2 Þ;

ðB:5Þ

@U ca c sinhðcb =2Þ2 ð1=2; 1=2Þ ¼  cba =caa b ; @x0 sinhðca Þ sinhðcb Þ

x¼0;1

Da ¼ DA =ðkA L2 Þ;

15

cbb ¼ cBB RB =kA :

Appendix B. Sufficient conditions for location and robustness of front

@W c sinhðca =2Þ2 cb ð1=2; 1=2Þ ¼  a þ cbb =cab ; @y0 sinhðca Þ sinhðcb Þ we find that if jUx ð1=2; 1=2Þj jcfa =caa j; jWy ð1=2; 1=2Þj jcfb =cab j, then jx0  1=2j; jy0  1=2j  1. Thus, for robustness to changes in f(x) for fixed values of ca ; cb ,

jcba =cfa j 1; jcbb =cfb j 1: Suppose there exists a solution ða1 ðxÞ; b1 ðxÞÞ to Eq. (A.2), with g ¼ H, the Heaviside function, and both a1 ; b1 are monotone decreasing. While this is not always the case, similar analyses can be performed otherwise. They are solved piecewise by

8   > x > > ; x 2 ½0; x0 ; < 1  k1 cosh pffiffiffiffi Da   a1 ðxÞ ¼ > > 1x > ffiffiffiffi ; x 2 ½x0 ; 1; : k2 cosh p

sinhðð1  y0 Þcb Þ ; sinhðcb Þ

With these parameter restrictions, we have the system

sinhðx0 ca Þ ; sinhðca Þ

k4 ¼

cab a1 ðy0 Þ þ cbb b1 ðy0 Þ þ rb þ cba f ðy0 Þ ¼ 0:

ðB:1Þ

Suppose a1 is symmetric about the midpoint when cfa ¼ 0; cfb ¼ 0, then x0 ¼ 1=2 and Eq. (B.1) becomes

1=2 þ cba =caa b1 ð1=2Þ þ r a =caa ¼ 0; a1 ðy0 Þ þ cbb =cab b1 ðy0 Þ þ r b =cab ¼ 0:

ðB:2Þ

For Eq. (B.2) to be invariant to changes in Db ; y0 must be equal to x0 . Therefore,

1=2 þ cba =caa 1=ð2lÞ þ r a =caa ¼ 0; 1=2 þ cbb =cab 1=ð2lÞ þ r b =cab ¼ 0; which gives

cba ¼ lð2r a þ caa Þ;

cbb ¼ lð2r b þ cab Þ:

ðB:3Þ

Define

Uðx0 ; y0 Þ ¼ aðx0 Þ þ cba =caa bðx0 Þ; Wðx0 ; y0 Þ ¼ aðy0 Þ þ cbb =cab bðy0 Þ: Then Eq. (B.1) is equivalent to

with c ¼ ð1 þ 2ra =caa Þ. Last, we will use the following change of ~ ¼ b  1=ð2lÞ, ~ ¼ a  1=2; b variables, a and define g~ðuA Þ ¼ gðcaa uA Þ  1=2 to get

ðB:7Þ

q ¼ cfa =caa ; h ¼ 1 þ 2 jcraaarj b . From here on, the tildes will be dropped

Db

caa a1 ðx0 Þ þ cba b1 ðx0 Þ þ r a þ cfa f ðx0 Þ ¼ 0;

at ¼ Da axx  a þ gðcaa ða  lcbÞ þ r a þ cfa f ðxÞÞ;   jra  rb a  lcb þ r b þ jcfa f ðxÞÞÞ bt ¼ Db bxx  lb þ gðjcaa 1 þ 2 caa j

~ þ qf ðxÞÞ; ~ þ g~ða ~  lcb ~ t ¼ Da a ~xx  a a ~ ~ ~ þ qf ðxÞÞÞ; ~ ~  lcb bt ¼ Db bxx  lb þ g~ðjðha

sinhðy0 cb Þ sinhðcb Þ

with ca ¼ p1ffiffiffiffi ; cb ¼ p1ffiffiffiffi. Therefore x0 ; y0 satisfy the system Da

ðB:6Þ

cab ¼ jð2r a þ caa Þ  2r b :

for constants x0 ; y0 ; k1 ; k2 ; k3 ; k4 . With continuity and differentiability of a; b enforced at x0 and y0 , we have

k3 ¼

cfb ¼ jcfa :

Using Eq. (B.3) with Eq. (B.6), we have

Db

k2 ¼

j :¼ cbb =cba ¼ cfb =cfa , then we have

cbb ¼ jcba ;

8   > > > l1 ð1  k3 cosh pxffiffiffiffiÞ ; x 2 ½0; y0 ; < Db   b1 ðxÞ ¼ > k > 1x > ffiffiffiffi ; x 2 ½y0 ; 1 : l4 cosh p

sinhðð1  x0 Þca ; sinhðca Þ

cba =cbb ¼ cfa =cfb : Define

Da

k1 ¼

Therefore it is natural to consider the two to be equal to a large factor, giving us the proportion

ðB:4Þ

in Eq. (B.7), viewing a; b as differences from a fixed concentration value and g as a shifted, scaled production function. Eq. (1) is obtained when h ¼ 1, or rb ¼ jra . References

[1] M. Bate, A. Arias, Development of Drosophila Melanogaster, Cold Spring Harbor Laboratory Press, 1993. [2] A. Ephrussi, D. St Johnston, Seeing is believing: the bicoid morphogen gradient matures, Cell 116 (2004) 143. [3] M. Klingler, The organization of the antero-posterior axis, Semin. Cell Biol. 1 (1990) 151. [4] L. Wolpert, Positional information and the spatial pattern of cellular differentiation, J. Theor. Biol. 25 (1) (1969) 1. [5] L. Wolpert, One hundred years of positional information, Trends Genet. 12 (1996) 359. [6] A.D. Lander, Morpheus unbound: reimagining the morphogen gradient, Cell 128 (2) (2007) 245. [7] W. Driever, C. Nusslein-Volhard, A gradient of bicoid protein in Drosophila embryos, Cell 54 (1988) 83. [8] W. Driever, C. Nusslein-Volhard, The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner, Cell 54 (1988) 95. [9] G. Struhl, K. Struhl, P.M. Macdonald, The gradient morphogen bicoid is a concentration-dependent transcriptional activator, Cell 57 (1989) 1259. [10] T. Berleth, M. Burri, G. Thoma, D. Bopp, S. Richstein, G. Frigerio, M. Noll, C. Nusslein-Volhard, The role of localization of bicoid RNA in organizing the anterior pattern of the Drosophila embryo, EMBO J. 7 (1988) 1749.

16

H. Hardway / Mathematical Biosciences 237 (2012) 1–16

[11] W. Driever, C. Nusslein-Volhard, The bicoid protein is a positive regulator of hunchback transcription in the early Drosophila embryo, Nature 337 (1989) 138. [12] R. Lehmann, C. Nusslein-Volhard, hunchback, a gene required for segmentation of an anterior and posterior region of the Drosophila embryo, Dev. Biol. 119 (1987) 402. [13] G. Struhl, P. Johnston, P.A. Lawrence, Control of Drosophila body pattern by the hunchback morphogen gradient, Cell 69 (1992) 237. [14] M. Simpson-Brose, J. Treisman, C. Desplan, Synergy between the hunchback and bicoid morphogens is required for anterior patterning in Drosophila, Cell 78 (1994) 855. [15] B. Houchmandzadeh, E. Wieschaus, S. Leibler, Establishment of developmental precision and proportions in the early Drosophila embryo, Nature 415 (2002) 798. [16] E. Lucchetta, J. Lee, L. Fu, N. Patel, R. Ismagilov, Dynamics of Drosophila embryonic patterning network perturbed in space and time using microfluidics, Nature 434 (2005) 1134. [17] T. Gregor, W. Bialek, R.R. de Ruyter van Steveninck, D.W. Tank, E.F. Wieschaus, Diffusion and scaling during early embryonic pattern formation, Proc. Natl. Acad. Sci. USA 102 (2005) 18403. [18] A.V. Spirov, D.M. Holloway, Making the body plan: precision in the genetic hierarchy of Drosophila embryo segmentation, Silico Biol. (Gedrukt) 3 (2003) 89. [19] D.M. Holloway, L.G. Harrison, D. Kosman, C.E. Vanario-Alonso, A.V. Spirov, Analysis of pattern precision shows that Drosophila segmentation develops substantial independence from gradients of maternal gene products, Dev. Dyn. 235 (11) (2006) 2949. [20] E.M. Lucchetta, M.E. Vincent, R.F. Ismagilov, A precise Bicoid gradient is nonessential during cycles 11–13 for precise patterning in the Drosophila blastoderm, PLoS One 3 (11) (2008) e3651. [21] Manu, S. Surkova, A.V. Spirov, V.V. Gursky, H. Janssens, A.-R. Kim, O. Radulescu, C.E. Vanario-Alonso, D.H. Sharp, M. Samsonova, J. Reinitz, Canalization of gene expression in the Drosophila blastoderm by dynamical attractors, PLoS Comput. Biol. 5 (3) (2009) e1000303. [22] Manu, S. Surkova, A.V. Spirov, V.V. Gursky, H. Janssens, A.-R. Kim, O. Radulescu, C.E. Vanario-Alonso, D.H. Sharp, M. Samsonova, J. Reinitz, Canalization of gene expression in the Drosophila blastoderm by gap gene cross regulation, PLoS Biol. 7 (3) (2009) e1000049. [23] T. Gregor, D. Tank, E. Wieschaus, W. Bialek, Probing the limits to positional information, Cell 130 (2007) 153. [24] O. Crauk, N. Dostatni, Bicoid determines sharp and precise target gene expression in the Drosophila embryo, Curr. Biol. 15 (2005) 1888. [25] F. He, Y. Wen, J. Deng, X. Lin, L.J. Lu, R. Jiao, J. Ma, Probing intrinsic properties of a robust morphogen gradient in Drosophila., Dev. Cell 15 (4) (2008) 558. [26] H. Hardway, B. Mukhopadhyay, T. Burke, T. James Hitchman, R. Forman, Modeling the precision and robustness of Hunchback border during Drosophila embryonic development, J. Theor. Biol. 254 (2008) 390. [27] J. Reinitz, D. Sharp, Mechanism of eve stripe formation, Mech. Dev. 49 (1995) 133. [28] L.F. Wessels, E.P. van Someren, M.J. Reinders, A comparison of genetic network models, Pac Symp. Biocomput. (2001) 508. [29] C. Schulz, D. Tautz, Autonomous concentration-dependent activation and repression of Kruppel by hunchback in the Drosophila embryo, Development 120 (1994) 3043. [30] T.J. Perkins, J. Jaeger, J. Reinitz, L. Glass, Reverse engineering the gap gene network of Drosophila melanogaster, PLoS Comput. Biol. 2 (2006) e51. [31] J. Treisman, C. Desplan, The products of the Drosophila gap genes hunchback and Kruppel bind to the hunchback promoters, Nature 341 (1989) 335. [32] J. Murray, Mathematical Biology, Springer-Verlag, 1989. [33] L. Edelstein-Keshet, Mathematical Models in Biology (Classics in Applied Mathematics), SIAM Society for Industrial and Applied Mathematics, 2005. [34] J. Jaeger, A. Martinez-Arias, Getting the measure of positional information, PLoS Biol. 7 (3) (2009) e81. [35] M. Howard, P. Rein ten Wolde, Finding the center reliably: robust patterns of developmental gene expression, Phys. Rev. Lett. 95 (20) (2005) 208103. [36] P. McHale, W.J. Rappel, H. Levine, Embryonic pattern scaling achieved by oppositely directed morphogen gradients, Phys. Biol. 3 (2006) 107. [37] B. Houchmandzadeh, E. Wieschaus, S. Leibler, Precise domain specification in the developing Drosophila embryo, Phys. Rev. E Stat. Nonlin. Soft. Matter Phys. 72 (2005) 061920. [38] A.M. Turing, The chemical basis of morphogenesis 1952, Bull. Math. Biol. 52 (1990) 153–197. [39] A. Gierer, H. Meinhardt, A theory of biological pattern formation, Kybernetik 12 (1972) 30. [40] H. Meinhardt, A. Gierer, Applications of a theory of biological pattern formation based on lateral inhibition, J. Cell. Sci. 15 (1974) 321.

[41] P. Gray, S.K. Scott, Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Isolas and other forms of multistability, Chem. Eng. Sci. 38 (1) (1983) 29. [42] J.E. Pearson, Complex patterns in a simple system, Science 261 (1993) 189. [43] A. Gierer, H. Meinhardt, A theory of biological pattern formation, Kybernetik 12 (1972) 30. [44] H.G. Othmer, E. Pate, Scale-invariance in reaction–diffusion models of spatial pattern formation, Proc. Natl. Acad. Sci. USA 77 (1980) 4180. [45] A. Hunding, P.G. Sørensen, Size adaptation of Turing prepatterns, J. Math. Biol. 26 (1988) 27. [46] T.C. Lacalli, L.G. Harrison, Turing’s conditions and the analysis of morphogenetic models, J. Theor. Biol. 76 (1979) 419. [47] B.N. Nagorcka, A pattern formation mechanism to control spatial organization in the embryo of Drosophila melanogaster, J. Theor. Biol. 132 (1988) 277. [48] A. Hunding, S.A. Kauffman, B.C. Goodwin, Drosophila segmentation: supercomputer simulation of prepattern hierarchy, J. Theor. Biol. 145 (1990) 369. [49] S.A. Kauffman, B.C. Goodwin, Spatial harmonics and pattern specification in early Drosophila development. Part II. The four colour wheels model, J. Theor. Biol. 144 (1990) 321. [50] B.C. Goodwin, S.A. Kauffman, Spatial harmonics and pattern specification in early Drosophila development. Part I. Bifurcation sequences and gene expression, J. Theor. Biol. 144 (1990) 303. [51] L. Diambra, L.d.a.F. Costa, Pattern formation in a gene network model with boundary shape dependence, Phys. Rev. E Stat. Nonlin. Soft. Matter Phys. 73 (2006) 031917. [52] T.C. Lacalli, D.A. Wilkinson, L.G. Harrison, Theoretical aspects of stripe formation in relation to Drosophila segmentation, Development 104 (1988) 105. [53] T.C. Lacalli, Modeling the Drosophila pair-rule pattern by reaction-diffusion: gap input and pattern control in a 4-morphogen system, J. Theor. Biol. 144 (1990) 171. [54] Y. Mori, A. Jilkine, L. Edelstein-Keshet, Wave-pinning and cell polarity from a bistable reaction–diffusion system, Biophys. J. 94 (2008) 3684. [55] Y.C. Wang, E.L. Ferguson, Spatial bistability of Dpp-receptor interactions during Drosophila dorsal-ventral patterning, Nature 434 (2005) 229. [56] D.M. Umulis, M. Serpe, M.B. O’Connor, H.G. Othmer, Robust, bistable patterning of the dorsal surface of the Drosophila embryo, Proc. Natl. Acad. Sci. USA 103 (2006) 11613. [57] F.J. Lopes, F.M. Vieira, D.M. Holloway, P.M. Bisch, A.V. Spirov, Spatial bistability generates hunchback expression sharpness in the Drosophila embryo, PLoS Comput. Biol. 4 (2008) e1000184. [58] M. Akam, Drosophila development: making stripes inelegantly, Nature 341 (1989) 282. [59] W.-M. Ni, P. Polacik, E. Yanagida, Monotonicity of stable solutions in shadow systems, Trans. Am. Math. Soc. 353 (2001) 5057. [60] Y. Nishiura, Global structure of bifurcating solutions of some reaction– diffusion systems, SIAM J. Math. Anal. 13 (4) (1982) 555. [61] J. Smoller, Shock Waves and Reaction–Diffusion Equations, Springer-Verlag, 1983. [62] A. Hillen, Classification of spikes and plateaus, SIAM Rev. 49 (1) (2007) 35. [63] T. Gregor, E. Wieschaus, A. McGregor, W. Bialek, D. Tank, Stability and nuclear dynamics of the bicoid morphogen gradient, Cell 130 (2007) 141. [64] T. Aegerter-Wilmsen, C. Aegerter, T. Bisseling, Model for the robust establishment of precise proportions in the early Drosophila embryo, J. Theor. Biol. 234 (2005) 13. [65] H. Jackle, D. Tautz, R. Schuh, E. Seifert, R. Lehmann, Cross-regulatory interactions among the gap genes of Drosophila, Nature 324 (1986) 668. [66] D.E. Clyde, M.S. Corado, X. Wu, A. Pare, D. Papatsenko, S. Small, A selforganizing system of repressor gradients establishes segmental complexity in Drosophila, Nature 426 (2003) 849. [67] J. Jaeger, The gap gene network, Cell. Mol. Life Sci. 68 (2011) 243. [68] M. Hoch, C. Schroder, E. Seifert, H. Jackle, cis-acting control elements for Kruppel expression in the Drosophila embryo, EMBO J. 9 (1990) 2587. [69] E. Segal, T. Raveh-Sadka, M. Schroeder, U. Unnerstall, U. Gaul, Predicting expression patterns from regulatory sequence in Drosophila segmentation, Nature 451 (2008) 535. [70] M. Hulskamp, C. Pfeifle, D. Tautz, A morphogenetic gradient of hunchback protein organizes the expression of the gap genes Kruppel and knirps in the early Drosophila embryo, Nature 346 (1990) 577. [71] R. Kraut, M. Levine, Mutually repressive interactions between the gap genes giant and Kruppel define middle body regions of the Drosophila embryo, Development 111 (1991) 611. [72] C.H. Waddington, Canalization of development and genetic assimilation of acquired characters, Nature 183 (1959) 1654.