| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
* Department of Applied Mathematics and Statistics and Center for Developmental Genetics, Stony Brook University, Stony Brook, New York 11794-3600
Department of Computational Biology, Center of Advanced Studies, St. Petersburg State Polytechnic University, St. Petersburg, 195251 Russia
Department of Biology, University of California, San Diego, California 92093
Universidade Federal do Rio de Janeiro, Instituto de Biofísica Carlos Chagas Filho, Rio de Janeiro, RJ 21949-900, Brazil
** Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545
1 Corresponding author: Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, NY 11794-3600.
E-mail: reinitz{at}odd.bio.sunysb.edu
 |
ABSTRACT |
|---|
 TOP ABSTRACT MATERIALS AND METHODSRESULTSDISCUSSIONACKNOWLEDGEMENTSLITERATURE CITED |
|---|
|
We illustrate the difficulties in interpreting mutant expression patterns with the example of the regulatory effect of Hb on Kr. The anterior boundary of the central Kr domain is shifted anteriorly in hb mutants (JäCKLE et al. 1986), while Kr expression is generally weaker than that in wild-type embryos (PANKRATZ et al. 1989). Moreover, embryos overexpressing hb show posterior expansion of the central Kr domain (HüLSKAMP et al. 1990). Finally, Kr expression is absent in embryos lacking both Bcd and Hb, but is restored in a concentration-dependent manner by reintroducing increasing dosages of Hb (STRUHL et al. 1992; SCHULZ and TAUTZ 1994). It has been proposed that these effects are due to a dual regulatory role of Hb with activation of Kr at low and repression of Kr at high concentrations of Hb (HüLSKAMP et al. 1990; STRUHL et al. 1992; SCHULZ and TAUTZ 1994).
However, the above observations can be explained equally well by indirect activation of Kr through Kni. The expression domain of kni, which encodes a repressor of Kr (JäCKLE et al. 1986; HOCH et al. 1992), expands anteriorly in hb mutants (HüLSKAMP et al. 1990), explaining reduced levels of Kr. The slightly altered posterior gt domain in hb mutants (ELDON and PIRROTTA 1991) further complicates interpretation, since Gt is a repressor of both Kr (KRAUT and LEVINE 1991b) and kni (ELDON and PIRROTTA 1991; CAPOVILLA et al. 1992). Expression of Kr is restored to high levels in hb kni double mutants (HARDING and LEVINE 1988), further supporting an indirect mechanism. Moreover, embryos overexpressing hb lack kni expression altogether (KRAUT and LEVINE 1991b), and posterior extension of the Kr domain in these embryos resembles Kr expression in kni mutants (JäCKLE et al. 1986). Finally, kni is widely expressed in embryos lacking Bcd and Hb, but is repressed in a concentration-dependent manner when Hb is reintroduced (STRUHL et al. 1992), which suggests that Kr derepression in these embryos is due to increasing repression of kni.
The above example reveals three main problems for inferring regulatory mechanisms from qualitative mutant expression data. These are the problems of consistency, uniqueness, and completeness.
Consistency of a proposed regulatory mechanism can be established only by keeping track of all regulatory inputs to a specific gene (cf. REINITZ and SHARP 1995). In the case of Kr, this involves at least five different regulatory contributions (by Bcd, Cad, Hb, Kni, and Gt). Current experimental approaches, however, are limited in their ability to monitor regulatory contributions simultaneously, as such interactions are inferred from mutant expression patterns and it is rarely possible to obtain mutants in more than three genes. Moreover, mutant regulatory systems by definition have an incomplete or otherwise defective set of regulatory interactions. Thus, the regulatory structure of the wild-type network must be assembled on the basis of evidence from different experiments. The consistency of such an inferred mechanism can be established conclusively only by testing it in the intact and complete developmental system.
Another problem for interpretation of mutant expression patterns is to establish the uniqueness of a mechanism, i.e., to decide whether regulatory interactions are direct or indirect. At least two possible regulatory mechanisms can account for the effect of Hb on Kr. Both mechanisms are consistent with available experimental evidence. In such an ambiguous situation, independent evidence can be provided by molecular approaches. Both Hb and Kni have been shown to bind to the Kr regulatory region in vitro (HOCH et al. 1991, 1992), but the functional importance of such biochemical interactions can be established only in vivo. Ideally, this would be achieved by targeted mutation of transcription factor binding sites in the regulatory region of an endogenous gene. Such an experiment is technically difficult and has not yet been attempted. Alternative approaches involving reporter constructs are subject to two significant complications. First, it is often difficult to establish the regulatory equivalence of such constructs to the endogenous gene. For instance, in kni mutants the posterior boundary of the third stripe of even-skipped (eve) is intact (FRASCH and LEVINE 1987), whereas the minimal enhancer for this stripe shows complete derepression between stripes three and seven (SMALL et al. 1996). Second, regulatory regions used in a construct may contain binding sites for multiple factors (see Kr above) or unknown binding sites, which leads to similar ambiguities in interpreting mutant expression patterns as in the case of the endogenous gene.
Finally, there is a fundamental issue concerning completeness of a proposed regulatory mechanism, which cannot be addressed by experimental approaches alone. The fact that all maternal and gap genes are necessary for correct gap gene expression does not prove that they are also sufficient. It is impossible to prove sufficiency of the inferred mechanism without reconstituting the system ab initio, using only well-defined ingredients. Such a reconstitution is obviously impossible by contemporary experimental methods and hence has to be attempted by using mathematical modeling and computer simulations.
The problems illustrated above show that to establish consistency, uniqueness, and completeness of a regulatory mechanism, we need a method that allows us to reconstitute wild-type gene expression patterns in silico, infer underlying regulatory interactions from these wild-type patterns, and keep track of all regulatory interactions in all nuclei at all times. The gene circuit method provides such an approach (MJOLSNESS et al. 1991; REINITZ and SHARP 1995; REINITZ et al. 1995, 1998). It is a data-driven mathematical modeling method whose main aim is to extract information about dynamical regulatory interactions between transcription factors from given gene expression patterns (Figure 2A; REINITZ and SHARP 1995). This is achieved in four steps: (1) formulation of a mathematical modeling framework, (2) collection of gene expression data, (3) fitting of the model to expression data to obtain regulatory parameters, and (4) biological analysis of the resulting gene circuits.
|
Gap gene circuits cover cleavage cycles 13 and 14A during the late syncytial blastoderm stage (Figure 2B; FOE and ALBERTS 1983), including most of embryonic stages four and five in CAMPOS-ORTEGA and HARTENSTEIN (1985). This covers the time between the first unambiguous detection of zygotically expressed Kr and Gt proteins in early cycle 13 (our own data and GAUL and JäCKLE 1987; ELDON and PIRROTTA 1991; KRAUT and LEVINE 1991a) and the onset of gastrulation at the end of cycle 14A (FOE and ALBERTS 1983). All nuclei divide equally and simultaneously at the beginning of cycle 14A.
Change in concentrations of transcription factors within each nucleus is governed by regulated protein synthesis, protein decay, and diffusion between neighboring nuclei (MJOLSNESS et al. 1991; REINITZ and SHARP 1995). Due to the lack of an in vitro polymerase II assay for eukaryotic transcription that faithfully reproduces in vivo transcriptional regulation, it is currently impossible to formulate a gene network model that is based on mechanistic chemical kinetics of transcription. Instead, we use coarse-grained kinetic equations for protein concentrations, which approximate the exact biochemistry with a sigmoid regulation-expression function (Figure 2C; MJOLSNESS et al. 1991; REINITZ and SHARP 1995).
Note that the general modeling framework outlined above does not specify which specific regulatory interactions take place within a gap gene circuit. These interactions are determined by regulatory parameters that constitute a genetic interconnectivity matrix (the T matrix). Each regulatory effect of a specific transcription factor on a target gene is described by a single parameter in the T matrix (Figure 2D). The gene circuit method aims to determine regulatory parameters and thus regulatory interactions within a gene circuit from given gene expression data. In other words, we seek sets of regulatory parameters that cause the gene circuit model to produce expression patterns that resemble real gap gene expression patterns as closely as possible (Figure 2A). This is achieved by fitting the model to quantitative segmentation gene expression data.
The set of quantitative gene expression data used in this study contains data for bcd, cad, hb, Kr, kni, gt, and tll from wild-type embryos (Figure 1; POUSTELNIKOVA et al. 2004). Data and model can be compared by numerically calculating expression patterns for given time classes from the model and then evaluating the sum of squared differences between model output and expression data for each gene, nucleus, and time class for which we have data. We minimize this sum by using a global optimization method called parallel Lam simulated annealing (PLSA, Figure 2A; CHU et al. 1999). The optimization procedure results in a gene circuit, which is defined by a specific set of regulatory parameters. Due to the stochastic nature of PLSA, different gene circuits (i.e., different sets of parameters) may be obtained, which all show essentially correct gene expression patterns.
The last step of the gene circuit method is the analysis of gene circuits to gain biological insights. The most important aspect of the gene circuit method considered here is that it allows for very detailed analysis of direct regulatory interactions within a given gene network. This is achieved by studying the distribution of gene circuit parameters between different gene circuits and by graphical analysis of regulatory contributions to specific patterning features (see RESULTS and REINITZ and SHARP 1995). This method of analysis allows us to study quantitative regulatory contributions to gene regulation in any nucleus at any point in time during a simulation.
Here we present a dynamical analysis of the gap gene network that is based on gap gene circuits. We show that these circuits are able to reproduce gap gene expression patterns in the late Drosophila blastoderm at high accuracy and temporal resolution. We provide a detailed analysis of regulatory interactions involved in gap gene regulation and show that our results are largely consistent with existing experimental evidence. Our models extend current knowledge of the gap gene system in several important aspects. We predict an activating effect of Cad on Kr and clarify evidence on the effects of Hb on Kr, Kr on kni, and Gt on kni. Our results suggest that mutual repression by complementary gap genes is absolutely essential for correct gap gene expression. We observe spatial asymmetry with posterior dominance in repressive interactions among overlapping gap genes. Moreover, the gene circuit method can provide information on regulatory mechanisms that is difficult to obtain by current experimental methods. Control of the posterior boundaries of posterior kni and gt was found to involve a temporal succession of multiple repressive interactions. Finally, we report a correlation between regulatory input from the terminal maternal system and late formation of gap gene domain boundaries in the posterior region of the embryo.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONACKNOWLEDGEMENTSLITERATURE CITED |
|---|
The internal state of nucleus i is described by concentrations vai of transcription factors encoded by segmentation genes denoted by index a. The change in transcription factor concentration over time, dvai/dt, depends on three processes during interphase: (1) protein synthesis, (2) protein diffusion, and (3) protein decay, represented by the summation terms on the right-hand side of Equation 1 below. During mitosis, protein synthesis is shut down and only diffusion and decay occur. Thus we write
 | (1) |
In Equation 1, Tab represents a matrix of regulatory coefficients where each coefficient Tab characterizes the regulatory effect of the product of gene b on the expression of gene a (Figure 2D). This matrix is independent of i, reflecting the fact that each nucleus contains a copy of the same genome. vBcdi is the concentration of Bcd in nucleus i. Bcd is exclusively maternal and its concentration is constant in time. The regulatory effect of Bcd on gene a is represented by the parameter ma. ha is a threshold parameter representing regulatory contributions of uniformly expressed maternal transcription factors. The relative rate of protein synthesis is then given by the sigmoid regulation-expression function 
, where 
is the total regulatory input on gene a (Figure 2C). The maximum synthesis rate for the product of gene a is given by Ra. The diffusion parameter Da(n) depends on the number of nuclear divisions n that have taken place before the current time t. Diffusion is assumed to vary inversely with the square of the distance between neighboring nuclei and this distance is halved upon nuclear division. 
a is the decay rate of the product of gene a. It is related to the protein half-life of the product of gene a by 
.
Quantitative expression data:
D. melanogaster blastoderm stage embryos were fluorescently stained for Eve protein and two other gene products using antibodies described in KOSMAN et al. (1998). As secondary antibodies, we used FITC anti-guinea pig, Texas Red anti-rabbit, and Cy5 anti-rat. Laterally oriented embryos were scanned using the 16x oil immersion objective of a Leica TCS4D confocal laser microscope. Fluorescent dyes were excited with a single wavelength at a time to ensure no leakage between channels, using the BP-FITC filter for the 488-nm excitation line (FITC), the BP-60030 filter for 568 nm (Texas Red), and the RG665 filter for 647 nm (Cy5). Expression levels were normalized per gene to a relative fluorescence intensity range of 0–255 on the basis of the most intensely fluorescent pattern on each slide with multiple embryos. Embryo images were cropped to fit embryo size and aligned along the A-P axis as shown in Figure 1.
Image segmentation:
A detailed description of this processing step can be found at http://flyex.ams.sunysb.edu/flyex/proc_steps/dave.html. Embryo images were segmented to obtain tabulated protein concentrations per nucleus as follows: Binary nuclear masks were constructed by edge detection, and average protein concentrations were obtained by averaging pixel values covering each nucleus in the mask. Nuclear positions are based on centroids of nuclei in the binary mask.
Time classification:
Embryos were assigned to cleavage cycle 12 (time class C12, used for initial conditions of the model at t = 0.0), cycle 13 (C13), and eight time classes (T1–T8) in cycle 14A (Figure 2B). Time classification for C12 and C13 is based on embryo morphology and for T1–T8 on careful visual inspection of the highly dynamic eve expression pattern by two independent observers (D. Kosman and S. Surkova; cf. MYASNIKOVA et al. 2001). Time classification was validated by membrane morphology (cf. MERRILL et al. 1988), as well as automated classification of eve expression patterns by complex-valued neural networks (AIZENBERG et al. 2002), and support-vector regression (MYASNIKOVA et al. 2002).
Background removal/registration:
Nonspecific background staining was approximated by a paraboloid and subsequently eliminated by a linear mapping of intensities that transforms fluorescence at or below background level to zero and transforms maximum fluorescence to itself (E. MYASNIKOVA, unpublished results). Expression patterns were registered using fast dyadic wavelets to align expression patterns as closely as possible (MYASNIKOVA et al. 2001). Only nuclei with positional values in the middle 10% along the D-V axis were used for further processing.
Integrated data:
Each integrated expression profile is based on registered data from at least 10 embryos stained for a specific gene at a specific time class, with the exception of Kni at C13, which is based on only two embryos, and Tll, for which we did not have data earlier than T3. Nuclei were categorized into 25 (C12), 50 (C13), and 100 (T1–T8) equal-sized bins according to their position along the A-P axis (cf. FOE and ALBERTS 1983). Concentration values for all nuclei in each bin were averaged to yield the final integrated one-dimensional expression pattern (Figure 1; POUSTELNIKOVA et al. 2004). The concentration of Bcd is nearly constant with respect to time during cycles 13 and 14A and is based on averaged registered bcd expression data from T1–T7. Concentrations of Cad and Hb at the onset of cycle 13 are derived from expression data for cycle 12. Initial concentrations for Kr, Kni, Gt, and Tll are zero in all nuclei.
Optimization by parallel Lam simulated annealing:
PLSA was used as described in REINITZ and SHARP (1995) and CHU et al. (1999). The set of ordinary differential Equations 1 was solved numerically using a Bulirsch-Stoer adaptive-step-size solver scheme adapted from PRESS et al. (1992). Equations were solved to a relative accuracy of 0.1%, and solutions were tested for numerical stability. We minimize the following cost function by adjusting parameters Ra, Tab, ma, ha, Da, and a in Equation 1:
 |
Summation is performed over the total number of data points Nd, i.e., the number of protein measurements across all genes a, nuclei i, and time classes t.
Parameter search spaces were defined by explicit search limits for Ra, Da, and a and a collective penalty function for Tab, ma, and ha as described in REINITZ and SHARP (1995). ha parameters of Kr, kni, gt, and hb were fixed to negative values representing a constitutive "off" state of the gene. This accelerated the annealing process considerably and slightly improved annealing results while not altering the overall quality of the resulting gene circuits. Optimization was performed in parallel on 10 2.4-Ghz Pentium P4 Xeon processors and took between 8 and 160 hr per optimization run.
Selection of gap gene circuits:
We use the root mean square (rms) score
 |
|
|
|
–0.005, (2) no interaction for parameter values between –0.005 and 0.005, or (3) activation for parameters 
0.005 (see Figure 4A). The threshold of 0.005 for the "no interaction" category was chosen empirically. Interactions falling into the no interaction category usually had no detectable effect on pattern formation in gap gene circuits analyzed graphically (see below). The gap gene network topology observed in a majority of gap gene circuits (Figure 4A) is preserved if a threshold of 0.01 is used instead (data not shown).
Software and bioinformatics:
Simulator and optimization codes were implemented in C; data quantification tools were implemented in C and the Khoros image analysis environment; and gene circuit analysis and plotting tools were implemented in Perl and Java. Software and gene circuit files are available at http://flyex.ams.sunysb.edu/lab/gaps.html. Expression data (FlyEx database) are available at http://urchin.spbcas.ru/flyex and http://flyex.ams.sunysb.edu/flyex.
|
RESULTS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONACKNOWLEDGEMENTSLITERATURE CITED |
|---|
|
|
Graphical analysis of gap gene regulation:
Graphical analysis of gap gene circuits allows us to "dissect" regulatory contributions of different transcription factors on the expression of a target gene and to characterize these interactions in great detail in space and time. To achieve this, we plot individual contributions to the sum of regulatory interactions affecting a gene's expression. Thereby, we focus on regions of expression domain boundaries. We identify regulatory factors responsible for the positioning of specific boundaries by looking for regulatory inputs that change significantly and consistently over the region of an expression domain boundary (cf. REINITZ and SHARP 1995). Consistent change implies that for boundary control by activation, the activator has to show a spatial expression gradient of the same polarity as the boundary it controls. Analogously, boundary control by repression implies a gradient of repressor with opposite polarity to the boundary it controls.
We have found activation of gap genes by Bcd and Cad in broad regions of the embryo (Figure 5). Bcd contributes strong activating inputs on the anterior domains of gt (Figure 5, A and C) and hb (Figure 5, E, F, H, and I) as well as on the central domain of Kr (Figure 5, B and D). Smaller activating inputs by Bcd can be detected in the posterior domains of kni (Figure 5, G and J) and gt (Figure 5, A and C). Three circuits (28003, 25005, and 29007) show repression of kni by Bcd, suggesting that Bcd activation might not be essential for kni expression during cycle 14A (Figure 4, A and C). The predominant maternal activating input on posterior kni and gt is provided by Cad (Figure 5, C and J). Furthermore, Cad provides a relatively strong activating input to central Kr expression (Figure 5D) and even contributes significantly to early anterior expression of hb (Figure 5H). Note that a small activating contribution of Cad on anterior hb can be detected in most gap gene circuits, but the strong early activation of hb by Cad shown in Figure 5H is exceptional. Activation in the posterior hb domain is largely due to Cad and hb autoactivation (Figure 4, A and E, and data not shown), a mechanism that we consider to be an artifact of the model (see DISCUSSION).
|
Whereas activation of gap genes by maternal genes occurs in rather broad regions, repressive interactions among gap genes provide spatially specific regulatory input for boundary positioning. Note that Kr and gt have mutually exclusive expression patterns in the blastoderm (Figure 1, G–I, and Figure 6A). Kr shows repression by Gt in all circuits (Figure 4, A and B). This repressive interaction is involved in positioning both anterior and posterior boundaries of central Kr expression during cycle 14A (Figure 6C). Although repression by Gt is quite strong, the regulatory profile of Kr indicates that missing repression by Gt does not lead to significant Kr derepression outside its central domain, since total regulatory input is not elevated significantly above the 10% level of expression in the absence of Gt (Figure 6C, arrow).
|
The anterior border of the posterior kni domain (Figure 1L and Figure 7, A and B) is set by a combination of repressive inputs by Hb and Kr (Figure 7, C and D). Whereas Hb represses kni in all circuits, repression by Kr was observed in only 6 of 10 circuits (Figure 4, A and C). Gap gene circuits without repression of kni by Kr show no detectable defects in kni expression (data not shown). Regulation of the posterior border of kni reveals a dynamic succession of repressive interactions during cycle 14A (Figure 7, E and F). All circuits show diminishing repressive input on kni by Tll in the region of the posterior boundary during cycle 14A (Figure 7, E and F), as tll expression is retained only in a region posterior to 80% A-P position (compare Figure 7A with 7B). In contrast, there is increasing repression by Gt and Hb in the boundary region (Figure 7, E and F). Note that the strength of repressive inputs by Gt and Tll varies greatly between circuits (Figure 4C and Figure 7, E and F). For instance, circuit 28008 (Figure 7E) shows extraordinarily strong repression of Gt, while other circuits such as 26001 show predominant repression by Hb and Tll, with a smaller contribution by Gt (Figure 7F).
|
|
|
DISCUSSION |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONACKNOWLEDGEMENTSLITERATURE CITED |
|---|
Some theoretical approaches to regulatory interactions in the segmentation gene network infer these interactions on the basis of interpretation of mutant expression patterns taken from the literature (SANCHEZ and THIEFFRY 2001; KUMAR et al. 2002). A difficulty with this approach is that such models tend to reproduce the interpretations of data they are based on, rather than providing an independent interpretation. In contrast, the gene circuit method does not require any a priori assumptions about specific regulatory interactions. Instead, it attempts to reconstruct these interactions on the basis of wild-type gene expression data (Figure 2A). Given this caveat, it is noteworthy that the results of our analysis of the gap gene network are largely consistent with studies based on mutant gene expression (see below). The fact that two independent methods lead to very similar results is an important cross-validation of conclusions based on both approaches.
Fitting models with many parameters to data is always at risk of producing nonspecific results. Gap gene circuits fail to fit expression data in regions of the embryo where additional factors are required for regulation, i.e., anterior of 
35% A-P position (cf. REINITZ et al. 1995), where gap gene regulation is known to involve head gap genes (COHEN and JüRGENS 1990; FINKELSTEIN and PERRIMON 1990; GROSSNIKLAUS et al. 1994), and posterior of 92% A-P position, where activity of the terminal gap gene hkb is required (WEIGEL et al. 1990; BRöNNER and JäCKLE 1991). Moreover, even though we have not obtained unique values for regulatory parameters in different circuits, we have found a strong tendency toward a specific type of regulatory interaction for most parameters (Figure 4). This suggests that gap gene circuits represent the gap gene regulatory network in a specific and reproducible way.
Mechanisms of gap gene regulation:
Although activating contributions from Bcd and Cad show some degree of localization (Figure 5), positioning of gap gene boundaries during cycle 14A is largely under the control of repressive gap-gap cross-regulatory interactions. Thereby, activation is a prerequisite for repressive boundary control, which counteracts broad activation of gap genes in a spatially specific manner (Figures 5–8). In addition, gap genes show a tendency toward autoactivation (Figure 4), which increasingly potentiates activation by Bcd and Cad during cycle 14A (Figure 5). Autoactivation is involved in maintenance of gap gene expression within given domains and sharpening of gap domain boundaries during cycle 14A. A similar, but less specific mechanism for spatially localized gene activation by maternal gradients has been proposed by MEINHARDT (1988).
Regulatory loops of mutual repression create positive regulatory feedback between complementary gap genes, providing a straightforward mechanism for their mutually exclusive expression patterns. Such a mechanism of "alternating cushions" of gap domains has been proposed by KRAUT and LEVINE (1991b) and CLYDE et al. (2003). Our results suggest that this mechanism is complemented by repression among overlapping gap genes. Overlap in expression patterns of two repressors imposes a limit on the strength of repressive interactions between them. Accordingly, repression between neighboring gap genes is generally weaker than that between complementary ones (Figure 4). Moreover, repression among overlapping gap genes is asymmetric, centered on the Kr domain (see Figure 9). Posterior of this domain, only posterior neighbors contribute functional repressive inputs to gap gene expression, while anterior neighbors do not. We show elsewhere that this asymmetry is responsible for anterior shifts of posterior gap gene domains during cycle 14A (JAEGER et al. 2004).
Repression by Tll mediates regulatory input to gap gene expression by the terminal maternal system (see Introduction). Tll provides the main repressive input to early regulation of the posterior boundary of posterior gt (Figure 8E), and activation by Tll is required for posterior hb expression (CASANOVA 1990; REINITZ and LEVINE 1990; BRöNNER and JäCKLE 1991). Note that these two features form only during cycle 13 and early cycle 14A (Figure 3), while other gap domain boundaries are already present at the transcript level during cycles 10–12 (KNIPPLE et al. 1985; TAUTZ 1988; MOHLER et al. 1989; ROTHE et al. 1989) and largely depend on the anterior and posterior maternal systems for their initial establishment (GAUL and JäCKLE 1987; TAUTZ 1988; MOHLER et al. 1989; RIVERA-POMAR et al. 1995). The delayed formation of posterior patterning features and their distinct mode of regulation are reminiscent of segment determination in primitive dipterans and intermediate germ-band insects, supporting a conserved dynamical mechanism across different insect taxa (TAUTZ and SOMMER 1995; DAVIS and PATEL 2002).
The set of regulatory interactions presented here provides a consistent and sufficient dynamical mechanism for gap gene expression (see Introduction). In summary, this set of interactions consists of the following five basic regulatory mechanisms (Figure 9): (1) broad activation by Bcd and/or Cad, (2) autoactivation, (3) strong repressive feedback between mutually exclusive gap genes, (4) asymmetric repression between overlapping gap genes, and (5) feed-forward repression of posterior domain boundaries by the terminal gap gene tll. In the following subsections, we discuss evidence concerning specific regulatory interactions involved in each of these basic mechanisms in some detail.
Activation by Bcd and Cad:
Activation of gap gene expression by Bcd and Cad is supported by the following. Bcd binds to the regulatory regions of hb, Kr, and kni (DRIEVER and NüSSLEIN-VOLHARD 1989; DRIEVER et al. 1989; HOCH et al. 1991; RIVERA-POMAR et al. 1995). The kni regulatory region also contains binding sites for Cad (RIVERA-POMAR et al. 1995). The anterior domains of gt and hb are absent in embryos from bcd mothers (TAUTZ 1988; ELDON and PIRROTTA 1991). The posterior domain of gt is missing in embryos mutant for both maternal and zygotic cad, while the posterior domain of kni is absent in embryos mutant for maternal bcd plus maternal and zygotic cad (RIVERA-POMAR et al. 1995). Our results suggest partial redundancy of activation of kni by Bcd, consistent with evidence from zygotic cad embryos from bcd mothers, where maternally provided Cad is sufficient to activate kni (RIVERA-POMAR et al. 1995).
Kr expression expands anteriorly in embryos from bcd mothers (GAUL and JäCKLE 1987), which is due to the absence of the anterior gt and hb domains (TAUTZ 1988; ELDON and PIRROTTA 1991; KRAUT and LEVINE 1991b). Bcd has been shown to activate expression of Kr reporter constructs (HOCH et al. 1990, 1991), supporting an activating effect of Bcd on endogenous Kr. The fact that Kr is still expressed in embryos from bcd mutant mothers has been attributed to activation by general transcription factors (KERRIGAN et al. 1991) or low levels of Hb (HüLSKAMP et al. 1990; STRUHL et al. 1992; SCHULZ and TAUTZ 1994). In contrast, our models predict that this activation is provided by Cad (Figure 4, A and B, and Figure 5D). Although Kr expression is normal in embryos overexpressing cad (MLODZIK et al. 1990), repressive control of Kr boundaries could account for the lack of expansion of the Kr domain in such embryos.
The activating effect of Cad on hb found in gap gene circuits is likely to be spurious. The anterior hb domain is absent in embryos from bcd mutant mothers (TAUTZ 1988), which show uniformly high levels of Cad (MLODZIK and GEHRING 1987). Moreover, the complete absence of the posterior hb domain in tll mutants (CASANOVA 1990; REINITZ and LEVINE 1990; BRöNNER and JäCKLE 1991) suggests activation of posterior hb by Tll rather than by Cad. We believe that this spurious activation of hb by Cad is due to the absence of hkb in gap gene circuits. The posterior hb domain fails to retract from the posterior pole in hkb mutants (CASANOVA 1990; BRöNNER and JäCKLE 1991), suggesting a repressive role of Hkb in regulation of the posterior hb border. Consistent with this, the posterior boundary of the posterior hb domain never fully forms in any of our circuits (Figure 3). Moreover, Tll is constrained to a very small or no interaction with hb (Figure 4E) due to the absence of the posterior repressor Hkb, since activation of hb by Tll would lead to increasing hb expression extending to the posterior pole.
Autoactivation:
A role for autoactivation in the late phase of hb regulation (SCHRöDER et al. 1988; HüLSKAMP et al. 1994) is supported by the fact that the posterior border of anterior hb is shifted anteriorly in a concentration-dependent manner in embryos with decreasing doses of zygotic Hb (SIMPSON-BROSE et al. 1994). Weakened and narrowed expression of Kr in mutants encoding a functionally defective Kr protein (WARRIOR and LEVINE 1990) suggests Kr autoactivation. Similarly, a delay in the expression of gt in mutants encoding a defective Gt protein (ELDON and PIRROTTA 1991) indicates gt autoactivation. However, our results suggest that gt autoactivation is not essential. It is generally weaker than autoactivation of other gap genes (Figure 4, B–E), and circuits lacking gt autoactivation show no specific defects in gt expression. Finally, in the case of kni, there is no experimental evidence for autoactivation, while some authors have even suggested kni autorepression (HOWARD 1990; ROTHE et al. 1994). We have not been able to detect such autorepression in any gap gene circuit (Figure 4, A and C).
Repression between complementary gap genes:
Mutual repression of gt and Kr is supported by the following. gt expression expands into the region of the central Kr domain in Kr embryos (ELDON and PIRROTTA 1991; KRAUT and LEVINE 1991a). In contrast, Kr expression is not altered in gt mutants before germ-band extension (GAUL and JäCKLE 1987; REINITZ and LEVINE 1990; ELDON and PIRROTTA 1991). However, Gt binds to the Kr regulatory region (CAPOVILLA et al. 1992), and the central domain of Kr is absent in embryos overexpressing gt (KRAUT and LEVINE 1991b). Moreover, Kr expression extends further anterior in hb gt double mutants than in hb mutants alone (KRAUT and LEVINE 1991b). The above is consistent with our analysis, which shows no significant derepression of Kr in the absence of Gt even though repression of Kr by Gt is quite strong (Figure 6C).
Hb binds to the kni regulatory region, and the posterior kni domain expands anteriorly in hb mutants (HüLSKAMP et al. 1990; ROTHE et al. 1994; CLYDE et al. 2003). Embryos overexpressing hb show no kni expression at all (NAUBER et al. 1988; ROTHE et al. 1989; KRAUT and LEVINE 1991b), and embryos misexpressing hb show spatially specific repression of kni expression (CLYDE et al. 2003). There is no clear posterior expansion of kni in hb mutants (HüLSKAMP et al. 1990; CLYDE et al. 2003). This could be due to the relatively weak and late repressive contribution of Hb on the posterior kni boundary or due to partial redundancy with repression by Gt and Tll (Figure 7, E and F). The posterior hb domain expands anteriorly in kni mutants, but anterior hb expression is not altered in these embryos (JäCKLE et al. 1986; CLYDE et al. 2003). Nevertheless, a role of Kni in positioning the anterior hb domain is suggested by the fact that misexpression of kni leads to spatially specific repression of both anterior and posterior hb domains (KOSMAN and SMALL 1997; WU et al. 2001; CLYDE et al. 2003). Moreover, only slight posterior expansion of anterior hb is observed in Kr mutants, while hb is completely derepressed between its anterior and posterior domains in Kr kni double mutants (CLYDE et al. 2003).
Repression between overlapping gap genes:
gt, kni, and Kr show repression by their immediate posterior neighbors hb, gt, and kni, respectively (Figure 4). Retraction of posterior Gt from the posterior pole during midcycle 14A fails to occur in hb mutants (MOHLER et al. 1989; ELDON and PIRROTTA 1991; KRAUT and LEVINE 1991a), and no gt expression is observed in embryos overexpressing hb (ELDON and PIRROTTA 1991; KRAUT and LEVINE 1991b). The posterior kni boundary is shifted posteriorly in gt mutant embryos (ELDON and PIRROTTA 1991), and kni expression is reduced in embryos overexpressing gt (CAPOVILLA et al. 1992). Note that these effects are very subtle and were not reported in similar studies by different authors (KRAUT and LEVINE 1991b; ROTHE et al. 1994). A weak but functional interaction of Gt with kni is consistent with our results. This interaction was found to be essential even in a circuit (29007) where it was deemed below significance level (Figure 4, A and C, and data not shown). Finally, Kni has been shown to bind to the Kr regulatory region (HOCH et al. 1992), and the central Kr domain expands posteriorly in kni mutants (JäCKLE et al. 1986; GAUL and JäCKLE 1987).
In contrast, we have been unable to detect any effect of Kr on hb (Figure 4, A and B). However, hb expression expands posteriorly in Kr mutants (JäCKLE et al. 1986; GAUL and JäCKLE 1989; CLYDE et al. 2003). This effect is likely to involve repression of hb by Kni. Kni levels are reduced in Kr embryos (PANKRATZ et al. 1989). hb is completely derepressed between its anterior and posterior domains in Kr kni double mutants, whereas anterior hb does not expand at all in kni mutants alone (CLYDE et al. 2003). Taken together with our results, this suggests that there is direct repression of hb by Kr in the embryo, but it is at least partially redundant with repression of hb by Kni.
Unlike repression by posterior neighbors, we have found no or only weak repression of posterior kni, gt, and hb by their anterior neighbors Kr, kni, and gt, respectively (Figure 4). Most gap gene circuits show weak activation of hb by Gt (Figure 4, A and E). Graphical analysis failed to reveal any functional role for such activation (Figure 5, H and I). Moreover, we have found no functional interaction between gt and Kni (Figure 4, A and D). Although relatively weak repression of kni by Kr was found in 6 out of 10 circuits (Figure 4, A and C), no specific patterning defects could be detected in the other 4. Consistent with the above, expression of posterior hb is normal in gt mutants, and both the anterior boundaries of posterior gt and kni are positioned correctly in kni and Kr mutant embryos, respectively (MOHLER et al. 1989; PANKRATZ et al. 1989; ELDON and PIRROTTA 1991; ROTHE et al. 1994).
Note that we have never observed activation of kni by Kr (Figure 4, A and C), which has been proposed to explain decreased expression levels of kni in Kr mutants (PANKRATZ et al. 1989; ROTHE et al. 1994). Our results strongly support the view that this interaction is indirect through Gt, which is further corroborated by the fact that kni expression is completely restored in Kr gt double mutants compared to that in Kr mutants alone (CAPOVILLA et al. 1992).
We have found a significant repressive effect of Hb on Kr (Figure 4, A and B). Consistent with this, Hb has been shown to bind to the Kr regulatory region (HOCH et al. 1991), and the central Kr domain expands anteriorly in hb mutants (JäCKLE et al. 1986; GAUL and JäCKLE 1987). However, partial redundancy of this interaction is suggested by correct positioning and shape of the anterior Kr domain in a circuit (28005) that does not show repression of Kr by Hb (Table 1).
It has been proposed that Hb plays a dual role as both activator and repressor of Kr (see Introduction). In the framework of the gene circuit model, concentration-dependent switching of regulative action could be implemented by allowing genetic interconnection parameters to switch sign at certain regulator concentration thresholds. Our current model explicitly does not include such a possibility. Nevertheless, we have been able to obtain circuits that reproduce Kr expression faithfully (Figure 3), suggesting that a dual role of Hb is not required for proper Kr expression. Moreover, we have never observed activation of Kr by Hb in any of the circuits (Figure 4, A and B). Therefore, our results support a mechanism in which the activation of Kr by Hb is indirect through derepression of kni.
Repression by Tll:
Only a few earlier theoretical approaches have considered terminal gap genes (MEINHARDT 1986; TCHURAEV and GALIMZYANOV 2001). Gap gene circuits accurately reproduce tll expression (data not shown). However, in gene circuits, tll is subject to regulation by other gap genes, which is inconsistent with experimental evidence (BRöNNER and JäCKLE 1991). In contrast, the correct expression pattern of tll in gap gene circuits allows us to study its effect on other gap genes in great detail. We have found strong repressive effects of Tll on Kr, kni, and gt (Figure 4). Tll binding sites have been found in the regulatory regions of Kr (HOCH et al. 1992) and kni (PANKRATZ et al. 1992). In tll mutants, Kr expression is normal (GAUL and JäCKLE 1987; REINITZ and LEVINE 1990), whereas expression of kni expands posteriorly (PANKRATZ et al. 1989), and the posterior gt domain fails to retract from the posterior pole (ELDON and PIRROTTA 1991; KRAUT and LEVINE 1991a). No expression of Kr, kni, or gt can be detected in embryos overexpressing tll under a heat-shock promoter (KRAUT and LEVINE 1991b; STEINGRIMSSON et al. 1991).
Comparison to logical analysis:
The logical analysis by SANCHEZ and THIEFFRY (2001) is the only other theoretical study of the gap gene system that achieves a level of detail comparable to the analysis presented here. Our results largely agree with SANCHEZ and THIEFFRY (2001), with the following notable exceptions. tll and the posterior hb domain were not considered in the logical analysis. The absence of posterior hb could explain why SANCHEZ and THIEFFRY (2001) did not report a r
