Multiple feedback loops achieve robust localization of wingless expression in Drosophila notum development

Journal of Theoretical Biology 292 (2012) 18–29

Contents lists available at SciVerse ScienceDirect

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

Multiple feedback loops achieve robust localization of wingless expression in Drosophila notum development Ken-ichi Hironaka a, Yoh Iwasa a, Yoshihiro Morishita a,b,n a b

Department of Biology, Faculty of Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan PRESTO, Japan Science and Technology Agency, 4-1-8 Honcho Kawaguchi, Saitama, Japan

a r t i c l e i n f o

a b s t r a c t

Article history: Received 15 February 2011 Received in revised form 1 August 2011 Accepted 20 September 2011 Available online 24 September 2011

Organ morphogenesis starts with the spatial patterning of different gene expressions in organ primordia, which is based on positional information provided by morphogens. To generate precise positional information, the robust localization of morphogen sources is needed. This can be realized by several different mechanisms, such as (i) by reducing the variations in spatial arrangement of morphogen sources, (ii) by reducing the variations in their source levels, and (iii) by increasing the degree of source localization with the sharp boundary that makes morphogen gradients steeper. Here we focus on the mechanism of localization of wingless expression, one of the important morphogens in Drosophila notum development. The mechanism of wingless-localization can be explained by a simple feed-forward loop network motif, but the real molecular network adopted by the organism is much more complex; it includes multiple feedback loops that function with the feed-forward loop in a coordinated manner. To clarify the functions of the molecular network, we decompose it into three submodules, each of which includes a single feedback loop, and examine their possible roles using a mathematical model. We demonstrate how the regulatory network for wingless expression realizes the conditions (i)–(iii) for its robust localization. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Systems biology Network motif Gene regulatory network Morphogen Positional information

1. Introduction Organ morphogenesis is performed in many steps in development. Typically it starts with the formation of spatial pattern of gene expressions in the organ primordia. In developing tissues, diffusive molecules called morphogens are secreted from their sources spatially localized in the tissues, and disperse through extracellular space by different ways such as diffusion and transcytosis. This generates spatial gradients of their concentrations (Entchev et al., 2000; Kerszberg and Wolpert, 1998; Lander et al., 2002; McDowell et al., 2001; Pfeiffer and Vincent, 1999). Each cell recognizes its own spatial position (i.e., the distance from its source) based on the concentrations of morphogens it receives, and regulates the expression levels of different genes according to the positional information (Wolpert et al., 2002). This step is followed by the construction of various anatomical structures, such as bones, muscles and bio-sensors (Gilbert, 1988; Wolpert et al., 2002), realizing tissue-specific morphogenesis. Generating precise positional information by morphogens requires proper spatial arrangement of morphogen sources as n Corresponding author at: Department of Biology, Faculty of Sciences, Kyushu University, 6-10-1 Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan. Tel.: þ81 92 642 2641. E-mail address: [email protected] (Y. Morishita).

0022-5193/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2011.09.022

well as an adequate level of morphogen production at each source (Fig. 1(a1) and (a2)). If the location of a morphogen source deviates from its proper position or if the amount of production is deviated from the adequate level, spatial profiles of the morphogen concentration change, leading to improper spatial patterning and to morphological anomalies. In addition, to give cells more accurate information of their position, a higher degree of source localization (or sharp boundary of the source) is desirable (Fig. 1(a3)). More localized source makes the morphogen gradient steeper. As reported in recent experiments with quantitative measurement, developmental processes would be affected by different kinds of noises (Blake et al., 2003; Elowitz et al., 2002; Gregor et al., 2007a, 2007b; Houchmandzadeh et al., 2002, 2005). Embryo–to–embryo variability of the source intensity (i.e., expression levels of morphogen molecules at the source) and embryo size can be the factors to prevent from generating precise positional information field in tissues. When morphogen gradients include such noises, corresponding deviation of the position is smaller for steeper gradient against the same magnitude of noises in morphogen concentrations (Eldar et al., 2003). If we quantify the precision of positional information in terms of the inverse of the difference between true position and position estimated by cells through morphogen concentrations, it is equal to the steepness of the morphogen gradient divided by the magnitude of noise in the

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

19

Fig. 1. (a) Generating precise positional information by morphogens requires reduced variations in spatial arrangement of morphogen sources (a1) and in the level of morphogen production at each source (a2). (a3) Spatially-concentrated source localization is also desirable for robust positioning. (b) Network motifs. (b1) (Incoherent type) feed-forward loop (FFL) in which the expression of common output genes (Z) are regulated by two or more signaling pathways (Y1 and Y2) whose activities are controlled by the same input signal (X). (b2) Positive (left) and negative (right) feedback loop (FBLs). (b3) A sub-network in which a FFL and a FBL are coupled. (c) Typical input–output curves of incoherent type FFL. (d) In the context of readout of morphogen gradient, an incoherent-type FFL can work as a detector of a specific distance from the source.

case of one-dimensional positioning by a single morphogen (Morishita and Iwasa, 2008, 2009). In the process of localization of morphogen sources, multiple morphogen species with different sources sometimes interact with each other. For example, in vertebrate limb bud development, fibroblast growth factor (FGF) with its source in the ectodermal tissue induces the localized expression of Shh within the mesenchymal tissue located with some distance from the FGF source. The induced Shh, in turn, activates Fgf expression (Laufer et al., 1994; Niswander et al., 1994). Another example is branching morphogenesis of the epithelial tube during lung development, in which SHH secreted from the tip of the tube induces the localized expression of Fgf in the mesenchyme, and while FGF ¨ 2006). signal activates Shh expression (Cardoso and Lu, In these reciprocal regulation between multiple morphogens, network motifs such as feed-forward loop (FFL) and feedback loop (FBL) are adopted in a coordinated manner to achieve robust

localization of morphogen sources against noises (Alon, 2007b; Milo et al., 2002). A FFL is a network motif in which the expression of common output genes are regulated by two or more signaling pathways whose activities are controlled by the same input signal (Mangan and Alon, 2003; Shen-Orr et al., 2002) (Fig. 1(b1)). A FFL is called ‘‘incoherent-type’’ when these signaling pathways have conflicting regulatory modes, i.e., one activating and the other repressing the output (Basu et al., 2004; Mangan et al., 2006). One important function of an incoherent-type FFL is to realize strong response only to the middle range intensity of input signal (Kaplan et al., 2008; Kim et al., 2008) (see Fig. 1(c)). In the context of responding to a morphogen spatial gradient, an incoherent-type FFL can work as a detector of a specific distance from the source (Fig. 1(d)). Such incoherent-type FFLs are often observed in many different systems including the spatial patterning of gap gene expressions in early development of Drosophila embryos, and the mesodermal patterning during gastrulation of

20

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

Xenopus blastula other than vertebrate limb and lung stated above. (Ishihara et al., 2005; Latinkic et al., 1997; Wolpert et al., 2002; Zeller et al., 2009) A feed-back loop (FBL) may be either negative or positive type (Fig. 1(b2)). Many experimental and theoretical studies on signal transduction systems have shown that a negative FBL has a noisecancelling function and a positive FBL can induce the bistability of the system, the latter realizing a switch-like response to the input signal of continuous intensity (Alon, 2007a; Laurent and Kellershohn, 1999; Savageau, 1974). Mechanisms for the localization of morphogen sources are likely to vary between organs and between species. For example, in the limb bud formation of vertebrates, one morphogen regulates the expression of another morphogen by the incoherenttype FFL, and the second morphogen in turn activates the first morphogen, which forms a positive FBL (Fig. 1(b3)). Hirashima et al. (2008) shows that a regulatory network of this type makes the distance between the two morphogen sources very robust against the fluctuation in the source intensity. In this paper, we mathematically study mechanisms for robust localization of wingless (wg) gene expression, one of the most important morphogen sources, in Drosophila notal development. As explained later, its expression is regulated by another morphogen, Dpp. In theory, the observed pattern of wg expression can be explained by a single FFL motif with Dpp signal as an input, but the real system adopted by the organism is much more complex. The main purpose of the present paper is to clarify possible functions of the complex system in which molecular components are studied well. We also discuss points necessary for experimental verification about whether or not the functions of the molecular network predicted by mathematical analyses work in the real system. The regulatory network for wg-expression is not specific for the notal development but adopted widely in other organs and in other species. Hence, the results of this study may contribute to the understanding of systems other than the Drosophila notal development. In the following, we first briefly review the molecular mechanism for Wingless (Wg) gene expression. Then we examine how robust localization of wg expression is achieved by combining FFL and FBL with mathematical models.

2. Molecular mechanism for the Dpp-dependent regulation of notal wg expression The notum is the dorsal portion of the mesothorax connected to wings (Fig. 2(a)). Dpp and Wg, which belong to the TGF-b and Wnt family, respectively, act as morphogens in the notal development as well as in the wing formation. Different structures including mechanosensory bristles and flight muscles develop in the notum as an outcome of cooperative action of Dpp and Wg (de Celis et al., 1999; Ghazi, 2003; Modolell and Campuzano, 1998). In the mutants, which have defects in Dpp signaling pathway, the expression patterns of wg and its downstream genes such as achaete–scute (ac–sc) and sr (stripe) (Ghazi, 2003; Tomoyasu et al., 1998) are very different from the normal patterns, resulting in the abnormalities in the muscle and mechano-sensor formations. In the whole wing disc, Dpp is expressed first along the A–P border which corresponds to the parasegment boundary in the embryonic ectoderm (Brook et al., 1996), and subsequently localized expression of wg becomes observed at a distance from the Dpp source (Fig. 2(b)). The notal wg expression is regulated by two signaling pathways in a Dpp-dependent manner (Garcia-Garcia et al., 1999; Sato and Saigo, 2000; Tomoyasu et al., 2000). One pathway represses wg expression through hetero-oligomer of transcription factors Pnr and Ush, whose expressions are induced by the Dpp

signal. The other pathway activates wg expression through Ushunbound Pnr. This regulation forms a typical incoherent-type FFL (Fig. 1(b1)). Therefore, the localized expression of wg at a distance from the Dpp source can be explained easily by the input–output relationship of the FFL. We can reproduce expression patterns of Pnr, Ush, and Wg observed in the experimental studies using a simple mathematical model for the FFL together with appropriate boundary conditions and reaction parameters (see Fig. 2(b) and Appendix A for the details of mathematical models). However, the regulatory network of wg expression adopted in Drosophila development is much more complex than this simple FFL (Fig. 2(c)). Multiple FBLs are imbedded in this network as explained below. In this paper, we attempt to clarify the role of these complex structures. We conjecture that coordinated operation of these FFL and FBLs might work to achieve robust localization of wg expression. We decompose the regulatory network of wg expression into three modules, each of which includes a single feedback loop. To understand the function of each feedback loop, we analyze the dynamics of each module in detail mathematically. To make clear the analysis and discussion, we here consider spatially-1D mathematical model along P–D axes (see Fig. 2(a)). In the main text, we explain only the key features of the model, but details of the model are given in Appendix B.

3. Mathematical analyses of three modules 3.1. Module 1: Mad phosphorylation by Dpp signal with negative feedback loop through Mad–Dad competition Module 1 is a commonly-used regulatory circuit in the Dpp signaling pathway, which appears in many developmental processes of Drosophila, such as establishment of embryonic dorsal–ventral polarity, gut formation, and patterning of imaginal discs (Massague´, 1998; Sekelsky et al., 1995) (Fig. 3(a)). The extracellular Dpp level is the input of this module, and it induces hetero-oligomerization of the type I receptor Tkv and the type II receptor Punt on the membrane. Punt activates Tkv by transphosphorylating the juxtamembrane region, and the activated Tkv phosphorylates intracellular signal mediator Mad (Smad family). The level of the phosphorylated Mad (pMad) is the output of the module. pMad forms homo-oligomeric complexes and/or hetero-oligomeric complexes with another Smad family protein Medea. Hetero-oligomers of pMad and Medea are transported to the nucleus where they transactivate many target genes. One such target, Dad, stably binds to Tkv and interrupts phosphorylation of Mad by Tkv, which forms a negative feedback loop in the Mad phosphorylation process (Inoue et al., 1998; Tsuneizumi et al., 1997). Biochemical reactions included in the module are expressed in terms of differential equations as Eq. (B1) in Appendix B1. Assuming that the time scale of those reactions is faster than that of tissue growth, we consider the input and output levels of the module in the steady state. After appropriate normalization (see Appendix B1 for details), we obtain the following input–output relationship: ½D ¼ eX ,   M* ¼

ð1aÞ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð½D þ a1 Þ2 þ 4f M ½D2 ð½D þ a1 Þ 2f M ½D

,

ð1bÞ

where [D] is the normalized Dpp level and [Mn] is the normalized pMad level (see Fig. 3(b) for the spatial profiles of [D] and [Mn]). X is the spatial coordinate normalized by the decay rate of Dpp spatial gradient. fM is the feedback strength, which is determined by quantities such as transcription rate of Dad by pMad and the dissociation constant of Dad–Tkv binding. When the transcription rate is larger and/or the dissociation constant is smaller, then the

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

Imaginal wing disc Notum Hinge Wing

21

Adult thorax Notum Hinge

Wing

Module1

Dpp Membrane Tkv-p

Tkv Mad

Mad-Tkv-p x

0

Proximal

Dad Dad-Tkv-p

Distal Mad-p Dad

b1

b2

b3

b4

Module2

[In vivo]

Ush

Pnr

b5

b6

b7

Pnr-Ush (complex)

b8

[In silico] Module3 Wgint Membrane Dpp (input)

Pnr

Ush

Wgext

Wg (output)

Fig. 2. (a) Imaginal wing disc (left), which is composed of prospective notum, hinge, and wing in adult thorax (right). (b) Comparison of spatial patterns of gene expressions in late third instar between in vivo and in silico. The left is proximal and the top is anterior. (b1–b4) Expression patterns of Dpp, Pnr, Ush, and Wg in vivo. pnr and ush is activated by Dpp signals secreted from the dorsal-most sources. wg expression domains are enclosed by white dots. Whereas pnr expression domain includes wg expression, ush expression domain is just next to wg expression domain. (b5)–(b8) Expression patterns in silico. Simple FFL model can reproduce similar spatial patterns in silico (see Appendix A for the details of mathematical models). (b1)–(b4) is minimally modified to serve this study, is reprinted from Mechanism of Development, vol. 93, pp. 127, Fig. 2A, C and E, Sato and Saigo, 2000 (with permission from Elsevier and from the corresponding author, Kaoru Saigo). (c) Schematic diagram for Dpp-dependent regulation of notal wg expression.

feedback strength becomes larger. a1 is the effective dissociation constant in the Mad phosphorylation. In the absence of the feedback loop (i.e., when fM is close to 0), we can simplify Eq. (1b) as follows:   M* ¼

½D : ½D þ a1

ð2Þ

In order to understand how the negative feedback loop through Mad–Dad competition affects the output variations against perturbations to the input, we calculate the sensitivity (Z) of the output to input variations, which is defined as Z  d ln½M*=d ln½D (Savageau, 1971). Smaller Z implies higher robustness against the perturbations. From Eqs. (1b) and (2), we can derive the following relation between the sensitivity in the presence of the feedback loop and the one in its absence:

a1 a1 ffio ¼ Znon-FB : ZFB ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð½D þ a1 Þ 2 2 ð½D þ a1 Þ þ 4f M ½D

effectively in the proximal region than in the distal region of the notum. We can conclude that Module 1 functions as a noisebuffer, which suppresses the propagation of the variations of the input to the downstream of signaling pathways.

ð3Þ

This indicates that the sensitivity in the presence of feedback loop is always smaller than that in the absence of feedback loop, regardless of parameters fM and a1. Further, inequality (3) clearly shows that the sensitivity decreases monotonically as the feedback strength fM increases (Fig. 3(c)). Fig. 3(d) shows the variations of the output at a certain position when the input signal includes variations following the Gaussian distribution. With the increase in fM, the variance of the output signal becomes smaller. The sensitivity is a monotonously decreasing function of the input level [D]. This suggests that the feedback works more

3.2. Module 2: Pnr and Ush expressions caused by pMad including a positive feedback loop via Pnr–Ush complex in Ush expression For Module 2, pMad works as the input signal and it activates its target genes, Pnr and Ush, as well as Dad. Pnr and Ush proteins form a complex, which functions as a repressor of wg expression. This complex also works as an activator of Ush expression, which forms a positive feedback loop in the Ush expression process (see Fig. 4(a1)). On the other hand, free Pnr that is not bound with Ush works as an activator of wg expression. We here examine the amounts of Pnr and Ush proteins, which are the outputs of this module. At the steady state, the input and output levels satisfy the following relationship (see Appendix B2 for details): ½P  ¼

½U  ¼

½M*h ½M*h þ ah2 ½M*h h

h 3

½M* þ a

,

ð4aÞ

þf U

ð½P½UÞh ð½P½UÞh þ ah4

,

ð4bÞ

22

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

Dpp

INPUT

Membrane 1 Tkv

fM = 0

Tkv-p

(without feedback)

Mad-Dad competition

Mad

Concentration

0.8

Dad Dad-Tkv-p

Mad-Tkv-p

fM = 1

0.6

fM = 10

0.4 Dpp

0.2 0 0

20

Mad-p

OUTPUT

40

60 80 100 Position (µm)

120

140

Dad

0.4 1

Frequency

Sensitivity

0.3

Dpp

0.6

fM = 10

Input

fM = 0

0.8

fM = 1

0.4

0.2

fM = 1

-0.2 -0.1 0 0.1 0.2

fM = 0

0.1 0.2

fM = 10 0

0 0

20

40

60

80

100

120

140

Position (µm)

-0.2

-0.1

0

0.1

0.2

Output variation

Fig. 3. Analysis of Module 1. (a) Schematic diagram for Mad phosphorylation process by Dpp signal with negative feedback loop through Mad–Dad competition. (b) Spatial profiles of p-Mad level with different feedback strengths (fM). (c) Spatial dependences of the sensitivity of p-Mad level against input variations for different feedback strengths. The sensitivity is smaller in more proximal regions of wing disc. The sensitivity becomes also smaller as the feedback strength increases. (d) Histograms of relative variations in pMad against input variations at x¼ 40 mm for different values of fM. Input variations were generated by Gaussian noise with the standard deviation s ¼ 0.3. Each histogram was made by 1000 numerical trials. With the increase in fM, the output variation becomes smaller, indicating that Module 1 functions as a noisebuffer which suppresses the propagation of the input variations to the downstream of signaling pathways. Parameter value: a1 ¼ 0.2 in all figures.

where [P] and [U] are the normalized levels of Pnr and Ush, respectively. In Eq. (4a) a2 gives the input level attaining halfmaximum of Pnr expression. By regarding a2 as a threshold of Pnr-expression, we determined its value from the position of distal boundary of Pnr expression observed experimentally. Since Pnr is an activator of wg expression, the distal boundary of the region for its expression also determines that of wg expression. Hill’s coefficient h represents the nonlinearity of the gene expression. In the following analyses, we mainly examine the case with h¼2, but we confirmed that the same conclusions hold for h ( 41) different from 2. In Eq. (4b) a3 and a4 are the input levels for half-maximum of Ush activation by pMad and by Pnr–Ush complex, respectively. In a similar way for Pnr expression, we determined the values of a3 and a4 using both positions of the distal boundary of Ush expression and the proximal boundary of wg expressions observed. fU indicates the strength of feedback by the complex in the Ush expression. An important role of a positive feedback loop in gene regulatory networks is to generate the bistability of systems, which make possible a switch-like (all-or-none) response of the system to continuously varying input levels (Ferrell Jr. and Xiong, 2001;

Hasty et al., 2001; Laurent and Kellershohn, 1999; Smolen et al., 2000). As shown in Fig. 4(b), the boundary of Ush expression is very sharp compared to the case without the feedback (i.e., the case of fU ¼0). Since Pnr–Ush complex is an repressor for wg expression, the sharp boundary of Ush expression makes the proximal boundary of wg expression sharp too. As stated earlier, localized expression of morphogen source with a sharp boundary is essential to provide precise positional information to cells in the developing field (Fig. 1(a3)). Furthermore, a sharp boundary of Ush contributes to forming a clear border between the medial and lateral notum, which are distinguishable by their cell adhesion properties (Calleja et al., 2000; Letizia et al., 2007). Hence the bistability appearing in this module is important in realizing robust normal development. To understand the system dynamics, we examined the bifurcation diagram of Ush expression level (i.e., an output of this module) for two parameters: the input level [Mn] and the feedback strength fU (see Fig. 4(a3) and (a4)). The parameter region with the bistability is wedge-shaped and is very narrow. Fig. 4(a2) shows input–output curves for different values of feedback strength. The bistability appears when fU is larger than a certain threshold value. For smaller input levels, i.e., to the left of the wedge-shaped bistable region in the diagram, the output stays small (in OFF state). In contrast, for larger

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

23

Fig. 4. (a) Analysis of Module 2. (a1) Schematic diagram for the regulations of pnr and ush expressions caused by pMad including a positive feedback loop of Ush expression via Pnr–Ush complex. (a2) Input–output curves of Module 2 for different feedback strengths (fU). Bistability occurs when fU is larger than a certain threshold. (a3) Bifurcation diagram for Module 2. A characteristic is a narrow and wedge-shaped bistable region, which makes possible switching response to input with high ON-OFF ratio (see the text for details). Parameter values: (a2,a3,a4) ¼(0.2, 0.4, 0.4). (a4) Parameter dependences of the bifurcation diagram on a3 and a4. (b) Spatial profiles of Pnr and Ush expression levels at the steady state. Ush profile becomes switch-like in the presence of the feedback loop. Parameter values: (a2,a3,a4,fU) ¼ (0.2, 0.4, 0.4, 3.0). (c1) Schematic diagram for regulations of pnr and ush expressions including a virtual positive feedback loop, auto-activation of Ush (Module 2n). This system was examined to compare the property of switching response between different feedback types. (c2) Input–output curves of Module 2n for different feedback strengths (fU0 ). (c3) Bifurcation diagram for Module 2n. In contrast to Module 2, the bistable region is wide. Especially, the system tends to have bistability even for small input. Parameter values: (a2, a3, a0 4) ¼(0.2, 0.4, 0.4). (c4) Parameter dependences of the bifurcation diagram on a3 and a0 4. (d) Comparison of sharpness of switching response between Module 2 and Module 2n, where the sharpness is defined as the difference between the maximum and the minimum of Ush expression levels in the bistable region (i.e., DU) divided by the width of input range for which the system has bistability (Dm).The right panel of Fig. 4(d) shows the histograms of the value of sharpness for the two types of FBLs for randomly chosen parameter sets (see Appendix C for details of calculation). Module 2 has a much higher sharpness than Module 2n.

input levels, the output becomes ON state. For middle levels of input, the output has two possible states, ON and OFF, both of which are locally stable. Which one should be realized depends on the initial

condition, history of the system, and noises experienced. This may cause the ambiguity of the boundary location of wg expression along the P–D axes, which is an undesirable situation for the purpose of

24

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

generating precise and robust positional information. We can guess that the narrow and wedge-shaped bistable region make possible to support the sharp switching between ON and OFF states against the input and less ambiguity of the boundary location of wg expression at the same time. To see if this property is specific to Module 2 with the feedback loop via the complex, we examined an alternative system (called Module 2n) with simpler mode of regulation—auto-activation of Ush as observed in many different gene regulatory networks (Fig. 4(c1)) (Milo et al., 2002; Shen-Orr et al., 2002; Thieffry et al., 1998). Fig. 4(c3) and (c4) shows the bifurcation diagram for the autoactivation system in which Eq. (4b) is replaced by the following one: ½U  ¼

½M* ½Uh 0 þf U , h h ½M* þ a3 ½U þ ða04 Þh

ð5Þ

where f0 U is the feedback strength and a0 4 is the level of Ush attaining the half-maximal effect of the feedback. In this auto-activation system (i.e., no Pnr–Ush complex involved), the parameter region for bistability is much broader than the system in Module 2 with positive feedback loop via Pnr–Ush complex. In Appendix C, we discuss mechanisms responsible for the difference in response between these two types of positive feedback loops. Here we calculate the sharpness of switching response for the two types, where the sharpness is defined as the difference between the maximum and the minimum of Ush expression level in the bistable region (see DU in Fig. 4(d)) divided by the width of input range for which the system has bistability (Dm). The right panel of Fig. 4(d) shows the histograms of the value of sharpness for the two types of FBLs for randomly chosen parameter sets (see Appendix C for details of calculation). Clearly Module 2 has a much higher sharpness than Module 2n. In this manner, we can conclude that the positive feedback loop via the Pnr–Ush complex realizes sharp switching between ON and OFF states against the input, which generates a clear (proximal) boundary of wg expression. Switch-like responses in biochemical reactions has been studied for the last 30 or 40 years and they are called ’’ultra-sensitivity’’ (Goldbeter and Koshland Jr., 1981; Koshland Jr., 1987). Typical systems with ultra-sensitivity are multiple-phosphorylation of a chemical and signaling cascade with multiple steps. The complexdependent feedback represented by the Module 2 studied in this section gives a new mechanism to realize ultra-sensitivity. 3.3. Module 3: wg expression regulated by Pnr and Pnr–Ush complex with auto-activation In Module 3, the input signals are Pnr and Pnr–Ush complex, and the output is the intracellular wg expression level or the extracellular concentration of its product (Fig. 5(a)). Pnr activates wg expression and the complex represses it. Since the boundary of Pnr expression is more distal than that of the complex, the bell-shaped expression of wg is achieved (Fig. 5(b)). The Wg product is secreted outside cells and diffuses through the extracellular environment, which provides positional information to cells around the source. It is reported that the extracellular Wg also up-regulates its own expression in the notum, although little is known about details of its molecular mechanisms (Sato and Saigo, 2000; Tomoyasu et al., 2000). This auto-activation is the third feedback loop in the system. In the steady state, the input and output levels satisfy the following: ! ah8 ½Ph ½W e h ½W i  ¼ 1þ f , ð6aÞ w h h h h ½P þ a6 ½W e  þ a7 ð½P½UÞh þ ah8 ½W e  ¼ a5

@2 ½W e  @X 2

þ ½W i ,

ð6bÞ

where [Wi] and [We] are the normalized levels of intracellular and extracellular Wg protein, respectively. In Eq. (6a) a6 is the input level for the half-maximal activation of wg expression by Pnr, and a8 is a similar quantity for the half-maximal repression by Pnr–Ush complex. These values can be estimated from the proximal and distal boundaries of wg expression observed experimentally. a7 is the input level for half-maximal effect in the auto-activation, and fw is the feedback strength. In Eq. (6b), a5 corresponds to the diffusion coefficient for the normalized extracellular Wg. As explained in the last section, the positive feedback loop promotes the switch-like response of the system. With the increase in fW, the region of wg expression is more localized, leading to a steeper gradient of the extracellular Wg protein level (see the top panel of Fig. 5(b)). In order to quantify the improvement of the precision of positional information given by Wg gradient realized by this feed-back loop, we calculated the spatial dependence of the precision, where the precision is defined as 9d[W]/dX9/s(X). The noise intensity s(X) was assumed to be proportional to its average value W(X). As shown in the bottom panel of Fig. 5(b), the precision is clearly improved as the feedback strength increases, especially on the proximal side of Wg expression region. This tendency can be more pronounced when we consider the case without the positive FBL in Module 2 (see Fig. 5(c)), suggesting that the functions of Module 2 and Module 3, the enhancement of localization of wg expression, may be redundant. We can conclude that Module 3 improves the localization of the expression of wg, which provides more precise positional information to the developing field.

4. Summary and discussion 4.1. Summary of the analyses To generate precise positional information field in tissues, robust localization of morphogen sources is an important process during organ morphogenesis. In this study, we focus on the localized gene expression of wg, one of important morphogens in Drosophila notum development, which is regulated by another morphogen Dpp. We decomposed the regulatory network for wg expression into three modules each of which includes a single FBL, and then analyzed the functions of each module one by one. Based on the analysis in this paper, we conclude that the robust localization of wg expression is performed through the following multiple steps. (i) First, wg expression occurs only in a region at a certain distance from the Dpp source due to the incoherent-type feed-forward regulation by the Dpp signal. (ii) Second, the Dpp signaling begins with the phosphorylation of Mad. The competitive inhibition of Mad phosphorylation by Dad forms a negative FBL (Module 1), and it works as noise-buffer that suppresses the propagation of variations of the input to downstream of signaling pathways. (iii) Third, Pnr and Ush, whose expressions are up-regulated by phosphorylated Mad, form a heterodimer and it positively regulates the Ush expression, which form a positive feedback loop (Module 2). This FBL make possible the switch-like response of Ush (or Pnr–Ush complex), which make the proximal boundary of wg expression steeper. (iv) Finally, auto-regulation of wg expression forms a positive FBL (Module 3), which also sharpens the boundary of wg expression. Consequently, steeper Wg gradient with proper source intensity and source location is realized, providing precise positional information to cells in a developing field.

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

25

INPUT Ush

Pnr

Pnr-Ush (complex)

INPUT

Auto-activation

Wgint

Membrane Wgext

OUTPUT

Pnr

1

Wgint

Dpp

0.8

Complex

0.6

Wgext (without FB) 0.4

Concentration

Concentration

1

Pnr

0.8

Dpp

0.6 0.4

0.2

0.2

0

0

0.08

fW = 10

fW = 5

0.08

0.06

Precision

Precision

Wgext

0.1

fW = 10

Wgext (without FB)

Complex

Wgext

0.1

Wgint

fW = 5

0.04

fW = 2

0.02

0.06 0.04 fW = 2

0.02

fW = 0

fW = 0

0

0 0

20

40

60

80

100

120

140

Position (µm)

0

20

40

60

80

100

120

140

Position (µm)

Fig. 5. Analysis of Module 3. (a) Schematic diagram for the regulation of wg expression by Pnr and Pnr–Ush complex with auto-activation. (b) (top) Spatial profiles of the amounts of Pnr, Pnr–Ush complex, intracellular Wg (with auto-activation of Wg), extracellular Wg (with auto-activation) and extracellular Wg (without auto-activation) at the steady state. The bell-shaped curve of wg expression can be explained by the input–output relationship of incoherent-type FFL. (bottom) Spatial pattern of precision of the positional information provided by Wg. The precision increases as the feedback strength fW increases. The increase of precision is more pronounced in proximal region. Parameters are(a5, a6, a8) ¼ (100, 0.3, 0.2),(fW,a7) ¼(0, any),(2.0, 0.08), (5.0, 0.2),(10.0, 0.44). For each value of fW, we adopted the value of a7 that gives the maximal precision of positional information. (c) Spatial profiles of the amounts of different chemicals (top) and spatial dependence of the precision of positional information in the case without the feedback in Module 2, i.e., fU ¼ 0 in Eq. (4b) (see the text). Parameter values: (a5, a6, a8)¼ (100, 0.3, 0.25). Parameter values for (fW,a7) are the same as in (b). The improvement of the precision by auto-activation of Wg is more pronounced in this case, suggesting that Module 2 and Module 3 are redundant with respect to their effect of enhancing the localization of wg expression.

4.2. Possible experimental verification on the function of each module conjectured by mathematical analyses Here we discuss possible ways for verifying the conjectured function of three modules. To verify the function of Module 1 as noise-buffer, it is necessary (i) to control chemical parameters that

change only the feedback strength fM keeping the effective dissociation constant a1 in the Mad phosphorylation, and (ii) to quantify the variations of pMad level at a focal location over different embryos (see Eqs. (1) and (3) and Fig. 3(d)). A candidate for such parameter is the transcription rate of Dad which are proportional to fM. The rate of de-phosphorylation of pMad is also a

26

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

candidate, which is inversely proportional to fM (see Appendix B1). Recently, it was reported that when the expression level of Tkv is altered experimentally, the variations of pMad in dad-deficient mutant become larger than those in the wild type (Ogiso et al., 2011). This result supports our results that Module 1 functions as a noise-buffer. We conjecture that the main function of Module 2 is to generate bistability, which enables a sharp boundary of Ush (and consequently Wg) expression. Especially, the characteristic of Pnr–Ush complex-dependent FBL is wedge-shaped bistability region in the bifurcation diagram (see Fig. 4a(3) and (4)), which is a clear difference from auto-activation system (see Fig. 4c(3) and (4)). Thus, to verify the function, making the bifurcation diagram experimentally by controlling feedback strength (fU) is necessary (see Eq. (4b)). Controlling the feedback strength is equivalent to changing the affinity between the complex and the promoter. At the present time, there is no report about possibility of controlling the affinity. Alternatively, the concept itself may be verified using synthetic systems such as Escherichia coli where a molecular circuit with the same topology as Module 2 is induced, because recent technology enables to construct artificially designed gene regulatory networks consisting of several genes (Andrianantoandro et al., 2006; Benner and Sismour, 2005). A bifurcation diagram similar to Fig. 4a(3) should support our prediction about the function of Module 2. To verify the function of Module 3 as sharpening the boundary of Wg expression, it is necessary (i) to block the auto-regulatory pathway in Wg expression and (ii) to quantify the amount of extracellular Wg proteins. As shown in the experiment by Sato and Saigo (2000), the auto-regulation is blocked in the mutant lacking dsh gene. Thus, by establishing the system to quantify the extracellular Wg, the verification should be possible. 4.3. Future extension of the model The model in the scheme of Fig. 2(c) does not include all the molecular elements (e.g., transcription factors) that are known to affect the notum development. For instance, it is reported that Bar homeobox interact with Dpp and Wg signal, and the role of this in spatial patterning needs to be clarified (Sato et al., 1999). In addition, if the Dpp signal controls the expression level of Dpp-receptor Tkv, which form another negative feedback loop, as observed in the development of wing pouch (Funakoshi et al., 2001), then the robustness of localization of wg expression in the notum may be improved further.

Acknowledgements This work was done by the support of a PRESTO procedure of JST to Y. Morishita, and a Grant-in-Aid for Scientific Research (B) of JSPS to Y. Iwasa.

Appendix A. Simple feed-forward loop (FFL) model can reproduce the expression patterns of Pnr, Ush, and Wg observed in experimental studies The FFL model we adopted in Fig. 2(b) is a spatially-2D model, which is given by the following partial differential equations: @q1 ¼ d1 r2 q1 g1 q1 , @t

ðA1aÞ

qh1 @q2 ¼ g2 q2 , h @t q1 þ K h21

ðA1bÞ

qh1 @q3 ¼ g3 q3 , h @t q1 þ K h31

ðA1cÞ

qh2 @q4 K h43 ¼ g4 q4 : h h h @t q2 þ K 42 q3 þ K h43

ðA1dÞ

Eq. (A1a) indicates the diffusion and degradation of input signal, Dpp, where its level is denoted by q1. d1 and g1 are the diffusion coefficient and the degradation rate. In the diffusion process of Dpp, we assume the constant source intensity whose value is set to 1 and the zero-flux boundary at the boundary except the source. Eqs. (A1b) and (A1c) indicate the dynamics of the activities of activation pathway (q2) and repressive pathway (q3) of wg expression (q4), both of which are regulated by the input. K21 (or K31) is the input level attaining the half-maximal activity of q2 (or q3), and g2 (or g3) is the degradation rate. Eq. (A1d) is for the output dynamics. K42 (or K43) is the level of q2 (or q3) attaining the half maximum in the activation (or repression) of wg expression. g4 is the degradation rate of Wg. The boundary shape of wing disc reflects the real one. We used the following parameter values in the simulation: d1 ¼1.0, g1 ¼0.002, g2 ¼1.0 , g3 ¼1.0, g4 ¼1.0, K21 ¼0.07, K31 ¼0.2, K42 ¼0.1, K43 ¼0.1.

Appendix B. Mathematical modeling of each module B1. Module 1 Biochemical reactions included in Module 1 is given as follows:     @½Tkv ¼ a1 ½Dpp Tkv þ b1 Tkv* , @t

ðB1aÞ

  @½Tkv*Mad ¼ a2 ½Tkv*½Madtot ða3 þ b2 Þ Tkv*Mad , @t

ðB1bÞ

    @½Mad* ¼ a3 Tkv*Mad b3 Mad* , @t

ðB1cÞ

     @½Tkv*Dad ¼ a4 Tkv* Dadtot b4 Tkv*Dad , @t

ðB1dÞ

  @½Dadtot  ¼ a5 Mad* b5 ½Dadtot  þc1 , @t

ðB1eÞ

½Tkv þ ½Tkv* þ ½Tkv*Mad þ ½Tkv*Dad ¼ T tot ðconstantÞ,

ðB1fÞ

Eq. (B1a) indicates the dynamics of the phosphorylation process of Tkv. The brackets [‘‘chemical name’’] indicate the concentrations of chemicals. Tkvn is the phosphorylated Tkv. The spatial gradient of Dpp is modeled by an exponential curve as ½Dpp ¼ D0 exp½lx, where D0 and l are the source level of Dpp and the decay rate of the gradient. x is the position along the proximal-distal axis, in which the proximal end is at x¼0. a1 and b1 are the phosphorylation and dephosphorylation rates, respectively. We assumed that total amount of Tkv (including free, phosphorylated and complex form) is constant (see Eq. (B1f)). Here, if the total amount of Tkv is larger than the total amounts of Mad and Dad, Tkvn Mad and Tkvn  Dad can coexist almost without competition. In order to verify the effect of Mad–Dad competition, we here assumed that the total amount of Tkv is much smaller than the total amounts of Mad and Dad, i.e., Mad and Dad in the free state are much larger than those in the modified states. Eq. (B1b) is for the dynamics of Tkvn  Mad complex. Madtot is the total amount of Mad. a2 and b2 are the association and dissociation rates, respectively. a3 is the phosphorylation rate. Eq. (B1c) is for the dynamics of the phosphorylation process. b3 is the dephosphorylation rate. Eq. (B1d) indicates the binding process of Dad with Tkvn. a4 and b4 are the association and dissociation rates,

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

27

respectively. Eq. (B1e) indicates the dynamics of the total amount of Dad. a5 is the transcription rate of Dad and b5 is the Dad degradation. c1 is the basal transcription rate regardless of the amount of phosphorylated Mad. In the steady state, the amount of normalized phosphorylated Mad, [Mn], satisfies the following equations,

In the steady state, the amounts of normalized Pnr and Ush proteins, [P] and [U], satisfy the following equations:

  M* ¼

½U  ¼

½D , ð1 þf M ½M*Þ½D þ a1

ðB1gÞ

  b2   M* ¼ Mad* , a a M T b1 ¼ 3 2 tot tot , b3 ðb2 þa3 Þ

ðB1iÞ

a2 M tot a4 c 1 þ , b2 þa3 b4 b5

ðB1jÞ

where [D] is the normalized Dpp level defined as [D]¼[Dpp]/D0. In the normalized system, we used the normalized spatial 2 coordinate X defined as X¼x/l. f M  a4 a5 b1 =b4 b5 b2 can be interpreted as the strength of negative feedback loop through Mad– Dad competition in binding with Tkvn, and a1  b1/a1b2 as the effective dissociation constant in the Mad phosphorylation as explained in the text.

Module 2 is modeled as follows: @½Pnr ½Mad*h ¼ a6 b6 ½Pnr @t ½Mad*h þ K h1 ða8 ½Pnr½Ushb8 ½Pnr2UshÞ,

ðB2aÞ !

  @½Ush ½Mad*h ½Pnr2Ushh ¼ a7 þf U b7 Ush h h h h @t ½Mad* þK 2 ½Pnr2Ush þ K 3

    @½Pnr2Ush ¼ a8 ½Pnr Ush b8 Pnr2Ush : @t

½M*h þ a2 h ½M*h ½M*h þ a3 h

,

ðB2dÞ

þf U

ð½P½UÞh ð½P½UÞh þ a4 h

,

ðB2eÞ

b6 ½Pnr, a6

ðB2fÞ

½U  ¼

 b7  Ush , a7

ðB2gÞ

where [Mn] is the normalized level of phosphorylated Mad. a2, a3, and a4 are defined as a2 ¼ b2K1/b1, a3 ¼ b2K2/b1, and a4 ¼b6b7b8K3/ a6a7a8, respectively (see the text for their meanings). B3. Module 3 Module 3 is modeled as follows: @½Wgint  ½Pnrh ¼ a9 @t ½Pnrh þ K h4

1 þf W

½Wgext h

½Wgext h þ K h5 b9 ½Wgint c4 ½Wgint  ,

!

K h6 ½Pnr2Ushh þK h6 ðB3aÞ

@½Wgext  ¼ DW r2 ½Wgext b10 ½Wgext  þ c4 ½Wgint : @t

B2. Module 2

ða8 ½Pnr½Ushb8 ½Pnr2UshÞ,

½M*h

½P  ¼ ðB1hÞ

b1

b2 ¼ 1þ

½P  ¼

ðB2bÞ ðB2cÞ

Eq. (B2a) indicates the dynamics of Pnr protein. a6 is the maximum rate of pnr expression induced by Madn, and K1 is the Madn level attaining the half-maximal activation in the pnr expression. b6 is the degradation rate. h is Hill’s coefficient representing the nonlinearity in the gene expression. We use the same symbol h to represent hill coefficients for different gene expression processes. We confirmed that main results in this study did not depend much on the magnitude of Hill coefficients, as far as they are greater than 1. The last two terms in the right hand sides are for the formation of the complex of Pnr and Ush. a8 and b8 are association and dissociation rates in the complex formation process. Eq. (B2b) describes the dynamics of Ush protein. The first term indicates the activation of ush expression by Madn. K2 is the Madn level attaining the half-maximal activation in the ush expression. a7 is the maximal rate of ush expression in the absence of the feedback through the Pnr–Ush complex. The second term indicates the activation of ush expression by the Pnr–Ush complex. fU is the feedback strength and K3 is the complex level attaining the half-maximal effect of the feedback. The feedback effect in the ush expression is modeled in a form of summation instead of product because non-zero expression of ush is observed when the binding between Pnr and Ush is inhibited. b7 in the third term is the degradation rate of Ush. Eq. (B2c) is the dynamics of the complex.

ðB3bÞ

Eq. (B3a) indicates the dynamics of the amount of intercellular Wg. a9 is the maximum rate of Wg expression in the absence of the feedback through Wg itself. K4 is the Pnr level attaining the half-maximal activation in the Wg expression by it. fW is the feedback strength and K5 is the level of extracellular Wg attaining the half-maximal effect of the feedback. K6 is the level of Pnr–Ush complex attaining the half-maximal repression in the Wg expression by it. b9 is the degradation rate of intracellular Wg and c4 is the transportation rate of the intracellular Wg to extracellular environment. Eq. (B3b) indicates the dynamics of the amount of extracellular Wg. DW is the diffusion coefficient of extracellular Wg, and b10 is its degradation rate. We adopted the zero-flux boundary condition in the Wg diffusion process. In the steady state, the normalized levels of intracellular and extracellular Wg products, [Wi] and [We], satisfy the following equations: ! ah8 ½Ph ½W e h ½W i  ¼ 1þf w , ðB3cÞ ½Ph þ ah6 ½W e h þ ah7 ð½P½UÞh þ ah8 ½W e  ¼ a5

@2 ½W e  @X 2

þ ½W i ,

ðB3dÞ

½W i  ¼

 b9 þ c4  Wgint , a9

ðB3eÞ

½W e  ¼

 b10 ðb9 þ c4 Þ  Wgext , c 4 a9

ðB3fÞ

where [P] and [U] are the amounts of normalized Pnr and Ush proteins. a5, a6, a7, and a8 are defined as a5 ¼DW/b10, a6 ¼b6K4/a6, a7 ¼b10(b9 þc4)K5/c4a9, and a8 ¼b6b7K6/a6a7, respectively (see the text for their meanings).

Appendix C. The difference in the diagrams between complexdependent FBL and auto-activation system As shown in bifurcation diagrams in Fig. 4(c3) and (c4), the bistable region for the system with auto-activation Eq. (5) is much

28

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

wider than that for the system with the complex-dependent feedback Eq. (4). This degrades the switching response of Ush expression to the input. In order to quantify the difference between the two types of feedback loops, we calculated the sharpness of switching response for the two types, where the sharpness is defined as the difference between the maximum and the minimum of Ush expression level in the bistable region (see DU in Fig. 4(d)) devided by the width of input range for which the system has bistability (Dm). The right panel of Fig. 4(d) shows the histograms of the value of sharpness for the two types of FBLs, where we randomly chose the set of parameter values for (a3, a4, fU) that realizes the left bifurcation point within [mn]¼ 0.170.01, and that are included in the ranges: a3A[0,1], a4A[0,1], and fUA[0,3]. As shown in the Figure, the complex-dependent feedback achieves higher sharpness than the auto-activation system. The difference can be understood as follows. When the input pMad level is non-zero, Eq. (4b) can be rewritten as follows: ½U  ¼

½M* ½Uh : þf 2 h h ½M* þ a3 ½U þ ðah4 =½Ph Þ

ðC1Þ

Eqs. (C1) and (5) differs in the threshold value at which the positive feedback loop starts to work. In the complex-dependent FBL, the threshold is a4/[P], a function of the input pMad level through the Pnr expression level. For Eq. (5), it is a constant a0 4. When the input is much smaller than a2 (see Eq. (4a)), the Pnr expression is very small and thus a4/[P] becomes very large, suggesting that FBL does not function at all. While, for larger input level, a4/[P] becomes closer to a0 4. This binary nature of the feedback effect makes possible to achieve the sharper switching response.

References Alon, U., 2007a. An introduction to systems biology: design principles of biological circuits. CRC Press. Alon, U., 2007b. Network motifs: theory and experimental approaches. Nat. Rev. Genet. 8, 450–461. Andrianantoandro, E., Basu, S., Karig, D.K., Weiss, R., 2006. Synthetic biology: new engineering rules for an emerging discipline. Mol. Syst. Biol., 2. Basu, S., Mehreja, R., Thiberge, S., Chen, M.-T., Weiss, R., 2004. Spatiotemporal control of gene expression with pulse-generating networks. Proc. Natl. Acad. Sci. U. S. A 101, 6355–6360. Benner, S.A., Sismour, A.M., 2005. Synthetic biology. Nat. Rev. Genet. 6, 533–543. Blake, W.J., Kaern, M., Cantor, C.R., Collins, J.J., 2003. Noise in eukaryotic gene expression. Nature 422, 633–637. Brook, W.J., Diaz-Benjumea, F.J., Cohen, S.M., 1996. Organizing spatial pattern in limb development. Annu. Rev. Cell Dev. Biol. 12, 161–180. Calleja, M., Herranz, H., Estella, C., Casal, J., Lawrence, P., Simpson, P., Morata, G., 2000. Generation of medial and lateral dorsal body domains by the pannier gene of Drosophila. Development 127, 3971–3980. ¨ J., 2006. Regulation of early lung morphogenesis: questions, facts Cardoso, W., Lu, and controversies. Development 133, 1611–1624. de Celis, J., Barrio, R., Kafatos, F., 1999. Regulation of the spalt/spalt-related gene complex and its function during sensory organ development in the Drosophila thorax. Development 126, 2653–2662. Eldar, A., Rosin, D., Shilo, B.-Z., Barkai, N., 2003. Self-enhanced ligand degradation underlies robustness of morphogen gradients. Dev. Cell 5, 635–646. Elowitz, M.B., Levine, A.J., Siggia, E.D., Swain, P.S., 2002. Stochastic gene expression in a single cell. Science 297, 1183–1186. Entchev, E., Schwabedissen, A., Gonz lez-Gait n, M., 2000. Gradient formation of the TGF-b homolog Dpp. Cell 103, 981–991. Ferrell Jr, J., Xiong, W., 2001. Bistability in cell signaling: how to make continuous processes discontinuous, and reversible processes irreversible. Chaos 11, 227–236. Funakoshi, Y., Minami, M., Tabata, T., 2001. mtv shapes the activity gradient of the Dpp morphogen through regulation of thick veins. Development 128, 67–74. Garcia-Garcia, M.J., Ramain, P., Simpson, P., Modolell, J., 1999. Different contributions of pannier and wingless to the patterning of the dorsal mesothorax of Drosophila. Development 126, 3523–3532. Ghazi, A., 2003. Prepattern genes and signaling molecules regulate stripe expression to specify Drosophila flight muscle attachment sites. Mech. Dev 120, 519–528. Gilbert, S., 1988. Development biology. Sinauer Associates. Goldbeter, A., Koshland Jr, D., 1981. An amplified sensitivity arising from covalent modification in biological systems. Proc. Natl. Acad. Sci. U. S. A 78, 6840–6844.

Gregor, T., Tank, D.W., Wieschaus, E.F., Bialek, W., 2007a. Probing the limits to positional information. Cell 130, 153–164. Gregor, T., Wieschaus, E.F., McGregor, A.P., Bialek, W., Tank, D.W., 2007b. Stability and nuclear dynamics of the bicoid morphogen gradient. Cell 130, 141–152. Hasty, J., McMillen, D., Isaacs, F., Collins, J.J., 2001. Computational studies of gene regulatory networks: in numero molecular biology. Nat. Rev. Genet. 2, 268–279. Hirashima, T., Iwasa, Y., Morishita, Y., 2008. Distance between AER and ZPA Is defined by feed-forward loop and is stabilized by their feedback loop in vertebrate limb bud. Bull. Math. Biol. 70, 438–459. Houchmandzadeh, B., Wieschaus, E., Leibler, S., 2002. Establishment of developmental precision and proportions in the early Drosophila embryo. Nature 415, 798–802. Houchmandzadeh, B., Wieschaus, E., Leibler, S., 2005. Precise domain specification in the developing Drosophila embryo. Phys. Rev. E 72, 061920. Inoue, H., Imamura, T., Ishidou, Y., Takase, M., Udagawa, Y., Oka, Y., Tsuneizumi, K., Tabata, T., Miyazono, K., Kawabata, M., 1998. Interplay of signal mediators of decapentaplegic (Dpp): molecular characterization of mothers against dpp, Medea, and daughters against dpp. Mol. Biol. Cell 9, 2145–2156. Ishihara, S., Fujimoto, K., Shibata, T., 2005. Cross talking of network motifs in gene regulation that generates temporal pulses and spatial stripes. Genes Cells 10, 1025–1038. Kaplan, S., Bren, A., Dekel, E., Alon, U., 2008. The incoherent feed-forward loop can generate non-monotonic input functions for genes. Mol. Syst. Biol. 4, 203. Kerszberg, M., Wolpert, L., 1998. Mechanisms for positional signalling by morphogen transport: a theoretical study. J. Theor. Biol. 191, 103–114. Kim, D., Kwon, Y.-K., Cho, K.-H., 2008. The biphasic behavior of incoherent feedforward loops in biomolecular regulatory networks. Bioessays 30, 1204–1211. Koshland Jr, D.E., 1987. Switches, thresholds and ultrasensitivity. Trends Biochem. Sci. 12, 225–229. Lander, A., Nie, Q., Wan, F., 2002. Do morphogen gradients arise by diffusion? Dev. Cell 2, 785–796. Latinkic, B.V., Umbhauer, M., Neal, K.A., Lerchner, W., Smith, J.C., Cunliffe, V., 1997. The Xenopus Brachyury promoter is activated by FGF and low concentrations ofactivinandsuppressed by high concentrationsof activin and by paired-type homeodomain proteins. Genes Dev. 11, 3265–3276. Laufer, E., Nelson, C.E., Johnson, R.L., Morgan, B.A., Tabin, C., 1994. Sonic hedgehog and Fgf-4 act through a signaling cascade and feedback loop to integrate growth and patterning of the developing limb bud. Cell 79, 993–1003. Laurent, M., Kellershohn, N., 1999. Multistability: a major means of differentiation and evolution in biological systems. Trends Biochem. Sci. 24, 418–422. Letizia, A., Barrio, R., Campuzano, S., 2007. Antagonistic and cooperative actions of the EGFR and Dpp pathways on the iroquois genes regulate Drosophila mesothorax specification and patterning. Development 134, 1337–1346. Mangan, S., Alon, U., 2003. Structure and function of the feed-forward loop network motif. Proc. Natl. Acad. Sci. U. S. A 100, 11980–11985. Mangan, S., Itzkovitz, S., Zaslaver, A., Alon, U., 2006. The Incoherent Feed-forward Loop Accelerates the Response-time of the gal System of Escherichia coli. J. Mol. Biol. 356, 1073–1081. Massague´, J., 1998. TGF-b SIGNAL TRANSDUCTION. Annu. Rev. Biochem. 67, 753–791. McDowell, N., Gurdon, J., Grainger, D., 2001. Formation of a functional morphogen gradient by a passive process in tissue from the early Xenopus embryo. Int. J. Dev. Biol. 45, 199–208. Milo, R., Shen-Orr, S., Itzkovitz, S., Kashtan, N., Chklovskii, D., Alon, U., 2002. Network Motifs: Simple Building Blocks of Complex Networks. Science 298, 824–827. Modolell, J., Campuzano, S., 1998. The achaete-scute complex as an integrating device. Int. J. Dev. Biol. 42, 275–282. Morishita, Y., Iwasa, Y., 2008. Optimal placement of multiple morphogen sources. Phys. Rev. E 77, 041909. Morishita, Y., Iwasa, Y., 2009. Accuracy of positional information provided by multiple morphogen gradients with correlated noise. Phys. Rev. E 79, 061905. Niswander, L., Jeffrey, S., Martin, G., Tickle, C., 1994. A positive feedback loop coordinates growth and patterning in the vertebrate limb. Nature 371, 609–612. Ogiso, Y., Tsuneizumi, K., Masuda, N., Sato, M., Tabata, T., 2011. Robustness of the Dpp morphogen activity gradient depends on negative feedback regulation by the inhibitory Smad,. Dad. Dev. Growth Differ. 53, 668–678. Pfeiffer, S., Vincent, J., 1999. Signalling at a distance: transport of wingless in the embryonic epidermis of Drosophila. Semin. Cell Dev. Biol. 10, 303–309. Sato, M., Saigo, K., 2000. Involvement of pannier and u-shaped in regulation of Decapentaplegic-dependent wingless expression in developing Drosophila notum. Mech. Dev. 93, 127–138. Sato, M., Kojima, T., Michiue, T., Saigo, K., 1999. Bar homeobox genes are latitudinal prepattern genes in the developing Drosophila notum whose expression is regulated by the concerted functions of decapentaplegic and wingless. Development 126, 1457–1466. Savageau, M., 1971. Parameter Sensitivity as a Criterion for Evaluating and Comparing the Performance of Biochemical Systems. Nature 229, 542–544. Savageau, M., 1974. Comparison of classical and autogenous systems of regulation in inducible operons. Nature 252, 546–549. Sekelsky, J.J., Newfeld, S.J., Raftery, L.A., Chartoff, E.H., Gelbart, W.M., 1995. Genetic Characterization and Cloning of Mothers against dpp, a Gene Required for decapentaplegic Function in Drosophila melanogaster. Genetics 139, 1347–1358. Shen-Orr, S.S., Milo, R., Mangan, S., Alon, U., 2002. Network motifs in the transcriptional regulation network of Escherichia coli. Nat. Genet. 31, 64–68. Smolen, P., Baxter, D., Byrne, J., 2000. Mathematical modeling of gene networks. Neuron 26, 567–580.

Ken-ichi Hironaka et al. / Journal of Theoretical Biology 292 (2012) 18–29

Thieffry, D., Huerta, A.M., Pe´rez-Rueda, E., Collado-Vides, J., 1998. From specific gene regulation to genomic networks: a global analysis of transcriptional regulation in Escherichia coli. Bioessays 20, 433–440. Tomoyasu, Y., Nakamura, M., Ueno, N., 1998. Role of Dpp signalling in prepattern formation of the dorsocentral mechanosensory organ in Drosophila melanogaster. Development 125, 4215–4224. Tomoyasu, Y., Ueno, N., Nakamura, M., 2000. The Decapentaplegic morphogen gradient regulates the notal wingless expression through induction of pannier and u-shaped in Drosophila. Mech. Dev. 96, 37–49.

29

Tsuneizumi, K., Nakayama, T., Kamoshida, Y., Kornberg, T., Christian, J., Tabata, T., 1997. Daughters against dpp modulates dpp organizing activity in Drosophila wing development. Nature 389, 627–631. Wolpert, L., Beddington, R., Jessell, T., Lawrence, P., Meyerowitz, E., Smith, J., 2002. Principles of development. Oxford University Press, Oxford. Zeller, R., Lopez-Rios, J., Zuniga, A., 2009. Vertebrate limb bud development: moving towards integrative analysis of organogenesis. Nat. Rev. Genet. 10, 845–858.