Reverse engineering the gap gene network of Drosophila melanogaster
A fundamental problem in functional genomics is to determine the structure and dynamics of genetic networks based on expression data. We describe a new strategy for solving this problem and apply it to recently published data on early Drosophila melanogaster development. Our method is orders of magnitude faster than current fitting methods and allows us to fit different types of rules for expressing regulatory relationships. Specifically, we use our approach to fit models using a smooth nonlinear formalism for modeling gene regulation (gene circuits) as well as models using logical rules based on activation and repression thresholds for transcription factors. Our technique also allows us to infer regulatory relationships de novo or to test network structures suggested by the literature. We fit a series of models to test several outstanding questions about gap gene regulation, including regulation of and by hunchback and the role of autoactivation. Based on our modeling results and validation against the experimental literature, we propose a revised network structure for the gap gene system. Interestingly, some relationships in standard textbook models of gap gene regulation appear to be unnecessary for or even inconsistent with the details of gap gene expression during wild-type development.
The segmented body pattern of Drosophila melanogaster is established by the expression of segmentation genes during the blastoderm stage of development (reviewed in [1,2]). The polarity of the embryo along its anterior–posterior (A–P) axis is established by maternal gradients of the Bicoid (Bcd), Hunchback (Hb), and Caudal (Cad) transcription factors (Figure 1A–1C). The trunk gap genes, hunchback (hb), Krüppel (Kr), knirps (kni), and giant (gt), are among the earliest targets of these maternal gradients. They show broad, overlapping expression domains (Figure 1E–1H). Gap genes together with maternal factors then regulate the expression of downstream targets, such as pair-rule and segment-polarity genes, which establish the segmental periodicity of the Drosophila body plan.
The gap gene network has been studied extensively by means such as the analysis of transcription factor binding sites on the DNA and the measurement of gene expression under wild-type and mutant conditions. These studies are the basis of qualitative regulatory models (e.g., [3,4]), which assert activating or repressing regulatory relationships between genes. These models provide a useful summary of the interactions within the network. However, ambiguities in the primary data can lead to different qualitative models. For example, if deleting a gene A results in decreased expression of a gene B, is A necessarily an activator of B? Or, could A be a repressor of some other gene C, which in turn represses B? This particular type of confusion has resulted in conflicting models of at least two relationships in the gap gene system, regarding the effect of Hb on Kr and the effect of Kr on kni. Qualitative models are also limited in that they usually do not specify precisely how conflicting regulatory signals (e.g., activation and repression) are resolved, leaving one without means for predicting the results of genetic or pharmacological interventions, for example. Further, qualitative models typically provide little or no explanation of the timing or dynamics of expression, which are key questions in processes such as development.
The increasing availability of quantitative gene expression data raises the possibility of detailed mathematical modeling of the regulatory relationships between genes. Indeed, recent work has shown that model parameters and even the existence of regulatory relationships can be automatically inferred from expression time series data [5–9]. Quantitative models can sometimes resolve ambiguities in qualitative models, offer much more specific predictions, and allow in-depth analysis of temporal changes in expression and the dynamical roles of regulatory relationships. For example, recent work on the gap gene system [6] has highlighted previously unappreciated effects of known regulatory relationships in shifting the expression domains of Kr, kni, and gt.
Unfortunately, fitting quantitative dynamical models can be computationally challenging. Only two techniques [5,10–12] have so far proven capable of fitting satisfactory models of the spatio–temporal segmentation gene expression data that we study, and both involve very long running times. Fitting all the models for the recent Jaeger et al. study [6,7], for example, which was done by parallel simulated annealing [5,10,11], required an estimated 20,000 CPU hours, or approximately two CPU years. Because of these intensive computing requirements, alternative regulatory explanations were not explicitly explored, and it has not been established whether previously proposed regulatory models (e.g., [3]) are equally capable of capturing the data. Long fitting times also make it difficult to explore alternative mathematical forms for the model. Gene regulation is an enormously complex process, and it is unclear in general how much of the biochemistry can be simplified away while retaining the key properties of the system. At a coarse level, the behavior of some genetic networks appears adequately captured by Boolean (on/off) models of gene expression [13–15]. However, detailed analysis of some genes has revealed more complex behavior which may well be functionally relevant [16,17]. Thus, it is advantageous to be able to easily experiment with different forms of models.
We use a novel strategy for fitting network models to spatio–temporal gene expression data that is compatible with a variety of modeling formalisms and, by combining function approximation techniques with simulation-based optimization, is vastly faster than previous methods [5,10–12]. The increased efficiency of the algorithm enables the fitting of many more models, allowing us to explore alternative modeling formalisms and assumptions about regulatory relationships between genes. We applied this strategy to fit four models to gap gene expression data [18]. Together, these models demonstrate different modeling formalisms and address alternative regulatory explanations for gap gene expression. Specifically, we show that both a smooth nonlinear representation of gene regulation (gene circuits [5–7]) and a discontinuous logical formalism based on activation and repression thresholds are capable of reproducing the observed gap gene expression dynamics, and that such models can efficiently be fit to the data. We show that it is possible to infer the activating and repressing relationships between genes directly from the data, as well as to test specific network structures, namely the model of Rivera-Pomar and Jäckle [3]. We find that the Rivera-Pomar and Jäckle model does not correctly capture all the major features of gap gene expression and requires additional regulatory links. Further, our optimizations consistently eliminate certain links in the Rivera-Pomar and Jäckle model, suggesting that those links are incorrect.
(In Protocol S1, we also report on fits using the more recent but similar Sanchez and Thieffry network structure [4], and find similar conclusions.) We summarize our findings in a Combined network model that includes a core of regulatory relationships found in all of our models in addition to regulatory relationships that are supported by some models and that are consistent with other experimental work.
Gap gene expression is established in the trunk region of the embryo, indicated by the black bars in Figure 1, during cleavage cycles 13 and 14A (before cellularization), a span of approximately 70 min. We model the wild-type expression of the four trunk gap genes, hb, Kr, kni, and gt (Figure 1C and 1E–1H). We do not model the expression of bcd, cad, or tailless (tll), a terminal gap gene (Figure 1A, 1B, and 1D). However, we allow the observed expression values for these genes to act as regulatory input to our models of the trunk gap genes. Expression of these genes is largely invariant along the dorsal–ventral axis of the embryo. Thus, we model expression as a function of time and of position along the A–P axis (Figure 1I–1L). At the start of cleavage cycle 13, there is some hb expression in the anterior, due to maternally deposited hb mRNA. There is no unambiguously detectable expression for the other three gap genes. By the end of cleavage cycle 14A, there are six main gap expression peaks, two each for hb and gt and one each for Kr and kni.
We represent gap gene dynamics using a reaction–diffusion partial differential equation [12]:
where va(x,t) denotes the relative concentration of gap protein a (unitless, ranging from 0 to 255) at space point x (from 35% to 92% of embryo length) and time t (0 min to 68 min after the start of cleavage cycle 13). The three terms on the right-hand side account for production, decay, and diffusion of protein, respectively. λa and Da are decay and diffusion rates. The model combines the processes of transcription and translation into a single production process. Pa specifies the production rate of protein a, as a function of the vector of concentrations of all proteins (including Bcd, Cad, and Tll) at the same point in space and time, v(x,t). The factor ζ(t) accounts for changes in transcribing gene density, due to the shutdown of transcription during mitosis [19] and due to the doubling of nuclei (see Materials and Methods for details).
We fit a total of four models, resulting from two different choices for the form of the production rate functions, Pa, and two different sets of constraints on the regulatory relationships between genes. The first model uses a gene circuit formalism identical to the one used by Jaeger et al. [6,7] for the production rate functions.
where Ra is the maximum production rate, b ranges over the seven genes {bcd, cad, hb, Kr, gt, kni, tll}, and . The regulatory weights, Tab, encode the effect protein b has on the production rate of protein a. If Tab > 0, then we interpret gene b as being an activator of gene a. If Tab < 0, then we say it is a repressor. We place no restrictions on the signs of the regulatory weights, so that the optimization determines the qualitative regulatory relationships between the genes. We call the resulting model Unc-GC, for unconstrained gene circuit. Following Jaeger et al. [6,7], we fix the bias, or offset, terms ha = −3.5 for all a.
While our first model tests the feasibility of inferring regulatory relationships de novo, our second model tests the ability of an established regulatory model to reproduce the data. In our second gene circuit model, regulatory relationships are limited to those in the model of Rivera-Pomar and Jäckle [3], with one exception. Rivera-Pomar and Jäckle did not consider the posterior hb domain, and their regulatory relationships cannot reproduce this domain (unpublished data). So, we added Tll activation of hb, which was sufficient to activate a posterior hb domain. We call this set of activation and repression relationships the RPJ network structure. The weights that represent these relationships are constrained to be positive or negative for relationships that are respectively activating or repressing. An additional complication in the RPJ network structure is that Kr is activated by low levels of Hb but repressed by high levels of Hb. We model this by allowing two weights for describing the effect of Hb on Kr—a positive weight that is multiplied by vhb(x,t) and a negative weight that is multiplied by (vhb(x,t))2/255. We call the resulting model RPJ-GC.
Previous analyses suggest that the general regulatory principle of the gap gene system is that genes are activated in broad regions by maternal gradients, but that repression from other gap genes can locally overwhelm general activation [3,7]. Our final two models rely on logical rules for the production rate functions that implement this general plan. Specifically, we assume that gap protein a is produced at rate Ra if at least one of its activators exceeds an activation threshold and none of its repressors exceed their repression thresholds. For example, the rule
states that Bcd is an activator of hb, hb is autoactivating, and Kr is a repressor. If there is enough Bcd or Hb and not too much Kr, as determined by the activation and repression thresholds, then there is production of Hb at the rate Rhb.
For our logical models, fitting the production rate functions means finding values for the Ra and for the activation and repression thresholds. We do not attempt to optimize which relationships are activating and which are repressing, although the optimization can change a threshold so that the activating or repressing effect is effectively eliminated. Our Unc-Logic model employs the regulatory structure discovered by our unconstrained gene circuit fit, except that we remove Gt activation of hb and Kni activation of gt (see Materials and Methods for justification). Our RPJ-Logic model uses the RPJ regulatory relationships.
We fit all four models using the same general strategy outlined in Figure 2 and described in greater detail in Materials and Methods and Protocol S2. In short, the strategy has three stages. Stage 1 produces an initial estimate of the decay and diffusion parameters for each gene and the spatial and temporal extents of protein production associated with each expression domain, without any regulatory explanation (Figure 2A and 2B). Stage 2 produces an initial estimate of the regulatory parameters by attempting to fit the estimated protein production rate at each space–time point (from Stage 1) in terms of the observed protein levels present (Figure 2C). In Stage 3, all model parameters are optimized using a local search procedure, starting from the initial regulatory, decay, and diffusion parameter estimates produced in the first two stages. The optimization approach we use in the third stage is much simpler than the approach used by Jaeger et al. [5–7,10,11]. However, if the initial parameter estimates are sufficiently good, then it can produce results of comparable quality and can do so relatively quickly.
The observed expression data is shown in Figure 3A and 3B. Simulated gap gene expression from our best-fitting models of each type is shown in Figure 3C–3F. The root mean squared errors (RMS error; see Materials and Methods) of these models are, in order of best to worst, 12.29 (Unc-GC), 14.83 (Unc-Logic), 15.88 (RPJ-GC), and 21.91 (RPJ-Logic). The Unc-GC and Unc-Logic models reproduce all six expression domains in approximately the right locations and with approximately the right timing (Figure 3C and 3D). The RPJ-Logic model suffers a major failure. It begins to form, but does not sustain, the posterior gt domain (Figure 3F). This accounts for much of the difference in RMS error between that model and the other four. The RPJ-GC model also contains a significant error, although it eventually reproduces all six domains. At the early time points, Kr production extends through the whole front of the embryo (Figure 3E, t = 11,24 min). Only later does repression from Hb eliminate production in the anterior, so that Kr is expressed in the correct region. The RPJ-GC model also creates a small erroneous Kr domain in the posterior half of the trunk (Figure 3E, t = 62 min). This is discussed further below.
The exact parameters for the best-fitting models of each type can be found in Protocol S3. The qualitative relationships they represent, as well as the relationships in some previously published models, are summarized in Figure 4A. Figures 5–8 provide detailed snapshots of the regulatory action between genes at different times. The models agree on many of the qualitative relationships between genes and never contradict each other, in the sense of one model claiming a relationship is activating while another claims it is repressing.
Figure 4 also includes a Combined network model and diagram, which represents our best estimate of the regulatory relationships in the real system. The Combined model includes a regulatory relationship if it was deemed significant in all of our models or if it was significant in it at least some of our models, performing a correct function, and is supported by other experimental evidence, such as binding site analysis or mutant or overexpression studies. Weights in a gene circuit model are considered insignificant if they are smaller than 0.0001 in magnitude. Terms in a logical model are considered insignificant if their removal results in no change in simulated expression. We adopted these stringent conditions on insignificance because we preferred to make decisions about weak/borderline interactions based on other experimental evidence rather than, for example, choosing some arbitrary threshold. (See the final paragraph of the Results and Discussion section for more information.) The regulatory structure of the Combined model is itself sufficient to reproduce all six gap gene domains using either the gene circuit or logical formalisms for production rate functions (Protocol S4). We first cite support for our Combined model, and then consider the results of the individual models in light of several outstanding questions about gap gene regulation.
The maternal proteins Bcd and Cad are largely responsible for activating the trunk gap genes, with Bcd being more important for the anterior domains and Cad more important for the posterior domains. Bcd is a primary activator of the anterior hb domain [20,21], the anterior gt domain [22], and the Kr domain [23,24]. Cad activates posterior gt [25]. The kni domain is present in bcd mutants and in cad mutants, but not in bcd;cad double mutants. This suggests redundant activation by the two maternal factors. Such redundant activation of kni is present in our Unc-GC model. For the other models, the optimization selected one or the other as activators, but not both. Tll is crucial for activating the posterior hb domain [26–28], while it represses Kr, kni, and gt, preventing their expression in the extreme posterior [22,29–34]. All the regulatory relationships between the gap genes in our Combined model are repressive. The complementary gap gene pairs, hb–kni and Kr–gt are strongly mutually repressive [22,32,35–38], as was found in nearly all of our models. (Repression of hb by Kni is not part of the RPJ regulatory relationships, but our Unc-GC and Unc-Logic models included the link.) Our models also suggest that mutual repression between hb and Kr helps to set the boundary between those two domains [38–41]. A chain of repressive relationships, hb⊣gt⊣kni⊣Kr [6,7], causes the shifts in the Kr, kni, and posterior gt domains. Autoactivation by hb is well-established [42], and there is also some evidence for autoactivation by Kr [43] and gt [22].
There is a long-running debate about whether or not low levels of Hb activate Kr. In hb mutants, the Kr domain expands anteriorly, suggesting that Hb represses Kr [39]. However, Kr expression in these mutants is lower than in wild-type [31] and expands posteriorly in embryos overexpressing Hb [36]. Further, in embryos lacking Bcd and Hb, the Kr domain is absent, but can be restored in a dosage-dependent manner by reintroducing Hb [44,45]. These observations suggest that Hb activates Kr. It has been suggested, therefore, that low levels of Hb activate Kr while high levels repress it [36,44,45]. An alternative explanation, however, is that the apparently activating effects of Hb are indirect, via Hb's repression of kni and Kni's repression of Kr [7]. Optimization of the Unc-GC model, which could have resulted in activation or repression of Kr by Hb, but not both, resulted in repression (Figure 4A). The RPJ models allow for a dual effect, but activation by Hb was eliminated during optimization of the RPJ-Logic model. The RPJ-GC model retained functional activation and repression of Kr by Hb. However, Kr expression in this model is defective. Kr is not properly repressed in the anterior (Figure 3E, t = 11,24 min). Further, Kr is ectopically expressed in a small domain in the posterior of the embryo (Figure 3E, t = 62 min and Figure 7C, t = 62 min). Thus, our models provide no support for activation of Kr by Hb. The only support we find, which is crucial in all models except Unc-Logic and also consistent with the mutant and overexpression studies cited above, is for repression of Kr by Hb.
Another point of disagreement in the literature is what prevents the expression of hb between its two domains. In the model of Rivera-Pomar and Jäckle [3], repression by Kr is the explanation. Our RPJ models confirm that this mechanism is sufficient. Specifically, in these models Kr repression prevents hb expression just to the posterior of the anterior hb domain (Figures 7B and 8B, t = 37,49,62 min). Between the Kr and posterior hb domains, there is no explicit repression of hb. Rather, Hb is not produced simply because of a lack of activating factors. In contrast, the models of Jaeger et al. [6,7] detected no effect of Kr and attributed repression solely to Kni. Our Unc-GC and Unc-Logic models found repression by Kni, but in addition to repression by Kr, not instead of it (Figures 5B and 6B). Kr is more responsible for repression near the anterior hb domain and Kni is more responsible for repression near the posterior hb domain. This is consistent with observations of expression in mutant embryos. Embryos mutant for Kr show slight expansion of the anterior hb domain [38], while kni embryos show expansion of the posterior hb domain [39]. In Kr;kni double mutants, hb is completely derepressed between its two usual domains [38]. This suggests, as seen in our Unc-GC and Unc-Logic models, that Kr and Kni are both repressors of hb, that their activity is redundant in the center of the trunk, and that Kr and Kni are the dominant repressors for setting the boundaries of the anterior and posterior domains, respectively. This interpretation was also favored by Jaeger et al. [6,7], on the basis of the mutant data, even though their models did not find repression by Kr.
In all of our models, the posterior hb domain is activated by Tll and sustained by Tll and hb autoactivation (Figures 5–8B, t = 37,49,62 min). Rivera-Pomar and Jäckle [3] did not consider the posterior hb domain, and did not include activation by Tll in their model. We added that link to the RPJ network structure because otherwise it was not possible to capture the posterior hb domain (unpublished data). The model of Jaeger et al. [6,7] captured the domain without Tll activation by substituting activation from cad. However, there is no confirming evidence for such an interaction. The absence of posterior hb in tll mutants [26–28] and the inability of our models to explain posterior hb by other means, leads us to believe the straightforward hypothesis that Tll activates posterior hb. Posterior hb is unique in that the domain begins to form later than the other five domains we model (Figure 3B, t = 37 min). In our RPJ models, this happens simply because high levels of Tll are needed to activate hb—levels that are reached only at about t = 30 min (Figure 3A, t = 24,37 min). The Unc-GC and Unc-Logic models also employ repression by Cad to slightly delay Hb production in the posterior (Figures 5B, t = 24 min, and 6B, t = 24 min). However, there is no confirming evidence for such repression, and we omit it from our Combined model.
Domain shifting was first observed by Jaeger et al. [6,7] and attributed to a chain of repressive regulatory relationships, hb⊣gt⊣kni⊣Kr. Our models largely support the importance of this regulatory chain, particularly the final two links. Repression of Kr by Kni was significant in all of our models (Figure 4A). Repression of kni by Gt was present in all models except RPJ-Logic, where it would be of little impact anyway, as RPJ-Logic has a defective posterior gt domain. Consistent with these findings, Kni binds to the regulatory region of Kr [29], and the Kr domain expands towards the posterior in kni mutants [39,40]. Similarly, the kni domain expands posteriorly in gt mutants [22], while embryos overexpressing gt show reduced kni expression [35].
Repression of gt by Hb is not as well supported by our models. The Unc-GC model included the link, though the regulatory weight was the smallest of all those in the model (Protocol S3). The link was eliminated from Unc-Logic and, of course, not present in the RPJ network structure. Instead, the models utilized decreasing activation by Cad (Unc-GC, Unc-Logic) and repression by Tll (Unc-GC, RPJ-GC) to shift the posterior gt domain (Figures 5–7D). Even with these links, however, shifting of the domain is not well-captured (Figure 3C, 3D, and 3F). RPJ-GC appears to capture the posterior gt shift best (Figure 3E). However, it relies on its small ectopic Kr domain to repress gt, a completely incorrect mechanism (Figure 7D, t = 62 min). Interestingly, a gene circuit fit using the network structure of Sanchez and Thieffry [4] (Protocol S1), captured the shift of posterior gt better than any of our other models, and it did so using repression of gt by Hb, providing additional modeling support for the relationship. There also is strong mutant evidence in favor of the relationship. In hb mutants, the posterior gt domain does not retract from the posterior pole [22,32,46]. Further, Gt is absent in embryos that have ubiquitous Hb, such as maternal oskar or nanos mutants [22,32] or embryos expressing Hb ubiquitously under a heat-shock promoter [33]. Thus, we find sufficient evidence to include a repressive link from hb to gt in our Combined model.
Activating or repressing links that oppose the direction of the repressive chain were eliminated by optimization of the Unc-Logic, RPJ-GC, and RPJ-Logic models (Figure 4A). In agreement with this result, the boundaries of the kni and posterior gt domains are correctly positioned in Kr and kni mutants, respectively [22,31,37,46]. Thus, the simplest picture supported by our models and consistent with the mutant studies is that there is no regulation from Kr, kni, or posterior gt to any of their immediate posterior neighbors, and that the repressive chain highlighted by Jaeger et al. [6,7] is indeed responsible for domain shifting.
All four of our models include autoactivation by hb. This is supported by the observation that late anterior hb expression is absent in embryos lacking maternal and early zygotic Hb [47]. Our models suggest hb autoactivation also plays a crucial role in sustaining the posterior domain, once it has been initiated by Tll (Figures 5–8B, t = 37–62), a role not previously emphasized. Autoactivation for the other genes was found by our Unc-GC model, but is not part of the RPJ network structure (Figure 4A). We included autoactivation only for Kr and gt in our Combined model, on the basis of a weakened and narrowed Kr domain in embryos producing defective Kr protein [43] and a delay in gt expression in embryos producing defective gt protein [22]. Interestingly, the gene circuit models of Jaeger et al. [6,7] also found autoactivation for all four gap genes, but they considered autoactivation by gt to be the weakest and least certain. In contrast, our Unc-Logic model retained gt autoactivation while eliminating autoactivation for Kr and kni (Figure 4A). The RPJ-Logic model was unable to reproduce the posterior gt domain. However, we found that by adding gt autoactivation to the model, it was able to create and sustain posterior gt correctly, bringing the error of the model down to 15.34 (unpublished data). This suggests that, after hb, gt is the most likely candidate for autoactivation. However, even this is not strictly necessary. The RPJ-GC model is able to reproduce and sustain the posterior gt domain without autoactivation by relying on cooperative activation from Bcd and Cad (Figure 7D).
The regulatory relationships proposed by Rivera-Pomar and Jäckle [3] are not fully consistent with the data and require amending. Repression of gt by Kni, which contradicts the mechanism of domain shifts described by Jaeger et al. [6,7], was eliminated by the optimization in both of our models based on the RPJ regulators. We also never observed activation of kni by Kr. We found no support for a dual regulatory effect of Hb on Kr. Activation of Kr at low levels of Hb was eliminated in the RPJ-Logic model. It was retained in the RPJ-GC model, but resulted in serious patterning defects. Inclusion of Tll as an activator of hb was sufficient to produce the posterior hb domain. Based on our fits and the primary experimental literature, there are likely other regulatory links missing from the model of Rivera-Pomar and Jäckle, though they are not strictly required to reproduce the wild-type gap gene patterns. Foremost is repression of hb by Kni, which appears important for eliminating hb expression anterior of the posterior domain. Fits based on the Sanchez and Thieffry regulatory relationships [4] (Protocol S1) also support these conclusions.
In contrast, the regulatory relationships in our Combined model and both the Unc-GC and Unc-Logic models are able to capture the wild-type gap patterns without gross defects. The relationships in the Unc-GC model are very similar to those obtained by Jaeger et al. [6,7]. For example, the regulation of Kr and kni is qualitatively equivalent in both models, and there is a single minor difference in the regulation of gt (Figure 4A). Our optimizations correctly identified activation of hb by Tll, which was missed by Jaeger et al. [6,7], though our models did less well at capturing shifting of the posterior gt domain. These regulatory relationships are also similar to those found by Gursky et al. [12], though that study was based on gap gene expression data with much lower accuracy and temporal resolution than the data used here [5,18]. These similarities show that differences in the mathematical formulations of these models—as ordinary versus partial differential equations, how diffusion and nuclei doubling are modeled, and choice of boundary conditions and other simulation parameters—are not important for the reproduction of the gap gene patterns nor for the inference of regulatory relationships from the data.