Reverse engineering the gap gene network of Drosophila melanogaster

Both of our gene circuit models fit the data better than the corresponding logical models. Although both types of models grossly simplify the complexity of gene regulation, this suggests that the gene circuit formalism is a better description of gap gene regulation than the logical formalism. However, the Unc-Logic model shows that the logical formalism can correctly capture the main features of gap gene expression. Of greater concern is that the strict on/off nature of the logical rules renders many regulatory inputs completely redundant, effectively eliminating them from the regulatory structure (Figures 4A, 6, and 8). In sparsely connected gene networks, this may be a useful bias. In a densely connected network like the gap gene system, it results in the elimination of many correct regulatory links. Another serious drawback of the logical models is that we could not find a satisfactory method for inferring the regulatory structure as part of the optimization process. We were forced to set the activators and repressors in the Unc-Logic model from the Unc-GC findings, and were able only to optimize the strengths of the regulatory relationships, represented by the activation and repression thresholds. Resolving this problem is an important avenue for further work.

The RMS errors in simulated expression for the Jaeger et al. [6,7] model and our Unc-GC model are comparable, so both optimization algorithms are equally successful in fitting the models to the data (Protocol S5). An important advantage of our technique, however, is speed. Each run of the algorithm used by Jaeger et al. [5–7,10,11] took approximately 3–10 days on a 10-processor machine, for a total of approximately two years of CPU time for their study. In contrast, it took on the order of a day or two to optimize each of our models, on a single-processor machine running unoptimized, uncompiled MATLAB code. The relative speed of our technique was crucial to our study because it allowed us to rapidly explore alternative modeling formalisms and to test specific network structures. The strength of our technique (see Materials and Methods) lies in the approach it uses to produce initial estimates of the regulatory, decay, and diffusion parameters. With good initial estimates, the models can easily be fined-tuned by a straightforward search algorithm. Indeed, the simplicity of the procedures we used to solve each of the three stages speaks to the power of the decomposition. However, it is certainly possible to substitute other, more sophisticated methods for solving each stage. For example, Stage 2 comprises a set of function approximation problems, for which many fitting techniques are available. The optimization approach of Gursky et al. [12] is quite fast when given a good starting point, and it would be interesting to use their algorithm for Stage 3 of our approach.

While our models, particularly Unc-GC, Unc-Logic, and the Combined model, capture the main features of gap gene expression dynamics, some failings are common to all the models. For example, none of the models capture well the shifting of the posterior gt domain (Figure 3). A gene circuit model based on the Sanchez and Thieffry regulatory relationships [4] (Protocol S1) does capture the shift, and by the expected mechanism, Hb repression of gt. However, the failure of the Unc-GC, Unc-Logic, and Combined models to capture the shift, despite including the necessary link, we can only attribute to imperfect optimization of the parameters. Our models also show defects in the early establishment of the domains, particularly the posterior gt domain (Figure 3C–3F, t = 11 min; similar defects were observed in [7]). This may be because the data is sampled less often early on, and so the RMS error effectively puts less weight on the correctness of the models early on. It may also be due to our simplification of transcription and translation into a single production process with no delay. At the start of cleavage cycle 13, gap gene transcripts begin to accumulate in the absence of gap proteins (except the Hb from maternally deposited mRNA). Thus, there should be no regulation between gap genes during the first few minutes of the time series. Our models allow regulation from the start of the time series, and this may be impacting the early expression patterns. Finally, none of the models capture the late parasegment 4-specific expression stripe of the anterior hb domain. This is visible as a small peak on top of the main anterior hb peak (Figure 1I or Figure 3B, t = 62 min). This expression stripe is due to a second hb promoter, different from the promoter responsible for the rest of the anterior hb domain [47].

It is likely that our gene circuit and logical models are simply incapable of capturing this phenomenon with a single production rate function, and would need to be generalized, perhaps by allowing hb production to be the sum of two separate production rate functions, to recreate the stripe.

Our models largely capture the wild-type expression patterns, and the regulatory relationships they rely on are consistent with mutant studies. However, our models do not in general display correct mutant expression patterns (see Protocol S6 for an example of simulating a Kr knockout mutant). The problem may point to a mismatch between biological reality and our modeling assumptions—for example, the mathematical forms of our production rate functions, our omission of production delays, or our treatment of protein concentrations as unitless numbers between 0 and 255. (Our expression data comes from images. Actual concentrations for the proteins are not known and may not even be of the same order of magnitude.) On the other hand, it is possible that the wild-type data to which the models were fit do not contain sufficient information to generalize to mutant organisms—although some gene circuit models fit to wild-type data do successfully predict mutant expression patterns [48]. Our fitting approach allows for rapid testing of alternative modeling assumptions, which may lead to models with improved predictions of mutant expression patterns. However, an important question for future research is whether we can find a way to incorporate the wealth of already available mutant expression data—which is often qualitative and not as temporally resolved as our wild-type data—into our fitting procedure.

Our models, particularly the Unc-GC model, includes a number of weak regulatory relationships, the significance of which is difficult to determine. Some of the links may arise from overfitting the data, or they may be compensating for incorrect modeling assumptions or for other missing or imperfectly modeled regulatory factors. However, it is difficult to say precisely which links should be ignored. Our data comprises a single space–time series, with strong correlations between datapoints. Resampling approaches for estimating significance, such as cross-validation or bootstrapping, are not useful in such cases. Jaeger et al. [7] deemed weights smaller in magnitude than a chosen threshold to be insignificant. However, the magnitudes of weights do not always reflect their importance to the model. For example, in the Unc-GC model, Tll activation of hb receives the smallest weight of any of hb's regulators, and yet this link is crucial for the formation of the posterior hb domain. In contrast, some links that receive larger weights in the Unc-GC model are dispensible, as shown for example by the success of the RPJ-GC model at reproducing the hb expression patterns. For these reasons, we used stringent criteria for dismissing links based on modeling alone. (Gene circuit weights of magnitude less than 0.0001 and logical terms that have absolutely no effect on expression were discarded.) Instead, we turned to external validation, admitting links to our Combined model only if they are additionally supported by other experimental evidence. However, it is also possible that some of the weak links in the Unc-GC or other models are correct. Regulatory effects between the gap genes can be subtle. For example, some studies have found changes in kni expression in gt underexpression and overexpression conditions [22,35], while other studies did not detect any change [33,37].

Additional, careful quantitative measurements of gap gene expression in wild-type as well as mutant organisms will be necessary to resolve the existence of the weaker links in our models.

Quantitative gene expression data used in this study are available online in the FlyEx database http://urchin.spbcas.ru/flyex or http://flyex.ams.sunysb.edu/flyex [18]. Methods for data acquisition and processing are described in detail elsewhere [49–52]. The data include quantitative wild-type concentration profiles for the protein products of bcd, cad, hb, Kr, kni, gt, and tll during cleavage cycles 13 and 14A, which constitute the late syncytial blastoderm stage of Drosophila development [53,54]. This covers a time of approximately 70 min between the first unambiguous detection of gap protein [7] and the onset of gastrulation [54]. Expression data from cleavage cycle 12 are used for the initial concentration of Hb at t = 0. Initial concentrations of Kr, Kni, and Gt are zero. The data represent expression at ten times, two during cleavage cycle 13 (t = 0 and 11 min) and eight during cleavage cycle 14A (t = 24,30,37,43,49,55,62, and 68 min). These times differ slightly from the ones used in [7], because we use a fixed time step of one minute for simulating our models (see below). The exact data are available as Dataset S1.

Here we describe our three-stage approach to fitting the parameters of the partial differential equations. Details of how this strategy was applied for each model can be found in Protocol S2. The ultimate goal is to produce a partial differential equation model [12]:

This is essentially the same as Equation 1, except that we have made explicit the regulatory parameters Θa that need to be optimized. Θa comprises Ra and either the weights of a gene circuit model or thresholds of a logical model. The term ζ(t) models the doubling of nuclei and shutdown of transcription during mitosis as follows.

From the experience of Jaeger et al. [6,7], it would appear that direct optimization of the parameters, Θa, λa, and Da is difficult. The first two stages of our approach are intended to produce good initial estimates of these parameters. The third stage is a direct optimization approach that fine-tunes the initial estimates to produce a good fit to the data (see Figure 2 for a summary).

In the first stage, we estimate λa, Da and the spatial and temporal extents of production associated with each domain. Peaks in gap protein expression are symmetrical, flat-topped, and have steeply sloping sides. It is reasonable to assume that, at any point in time, the rate of protein production associated with a particular peak is constant in some interval along the A–P axis and zero outside of that interval. However, it is not clear from the data exactly where such boundaries should be. It is also not clear exactly when production begins and when (or if) it ends.

We use a set of seven parameters to describe the conditions of protein production associated with each of the six gap protein peaks: ρ, a production rate; τstart, the time at which production begins; τend, the time at which production ends; xs,a, the anterior-most extent of production at time τstart; xs,p, the posterior-most extent of production at time τstart; xe,a, the anterior-most extent of production at time τend; and xe,p, the posterior-most extent of production at time τend. The last six parameters define a quadrilateral region of space–time (Figure 2A). Within this quadrilateral, we take protein production to occur at rate ρ. Genes with multiple peaks are assigned one quadrilateral per peak. At space–time points that do not fall in any of the quadrilaterals for a particular protein, the rate of protein production is assumed to be zero. Thus, for a single-domain gene, the production rate function has the form

For a two-domain gene, the production rate function is the maximum of two such functions. In either case, the dynamical model of the expression of protein a is

where are parameters that specify the spatial extent of production and production rates for the expression domain or domains of protein a. To fully specify Equation 7 requires nine parameters for a single-domain gene and 16 parameters for a two-domain gene—a decay rate, a diffusion rate, and seven parameters describing the production conditions for each domain. The system of partial differential equations is decoupled—the genes do not interact—so the parameters for each gene can be optimized independently. We optimize the parameters for each gene to minimize the RMS error between simulated and observed expression, , where va(x,t) is the simulated expression of protein a, ya(x,t), is the observed expression, and Nd is the number of space–time points. We used a standard, repeated first-improvement local search with randomized order of neighbor examination to optimize the parameters for each gene. Stage 1 is performed only once for all four models, as it does not depend on the form of the Pa or the regulatory structure or parameters.

Next, we generate an initial estimate of the regulatory parameters for each gene by searching for Θa that minimize the error function , where Nd is the number of space–time points, is the Stage 1 estimate for the production rate of protein a at space–time point (x,t), and Pa is the production function for protein a (as in Equation 4, but with the observed expression values given as input). In other words, we search for regulatory parameters which fit, as closely as possible, the quadrilateral production regions found in Stage 1. Parameters can be optimized separately for each gene. This is just a problem in function approximation (also known as regression or supervised learning)—we seek a set of parameters Θa so that Pa reproduces as well as possible the input–output pairs . Many techniques are available for solving this sort of problem. For our gene circuit models, the error function is differentiable with respect to Θa. We optimize the parameters using repeated runs of an adaptive step-size gradient descent procedure. For our logical models, the regulatory parameters are optimized using a repeated first-improvement local search with randomized order of neighborhood examination.

Although the Unc-GC model includes Gt activation of hb and Kni activation of gt, we disallowed these links in the Unc-Logic model. There is no support in the literature for these links and we found that their removal improved our fits.

Finally, we combine the decay and diffusion constants estimated in Stage 1 with the regulatory parameters estimated in Stage 2 in a fully coupled partial differential equation model (Equation 4). Starting from these initial parameters, we perform repeated first-improvement local search with randomized order of neighbor examination, seeking parameters that minimize the RMS error , where Nd is the number of genes modeled (4) times the number of space–time points, va(x,t) is the simulated expression from the model, and ya(x,t) is the observed expression.

We simulate using a fixed time step of one minute and a spatial grid of 58 points (one space point for each 1% of embryo length between 35% and 92%). For each step, we calculate the production rates for each gene at each space point and add them to the expression values. This corresponds to one minute of constant-rate production with no decay or diffusion. For Equation 1 (or 4), in which production rates depend on the protein levels present, the simulated gap protein levels and the observed bcd, cad, and tll levels are used to calculate production. Values of bcd, cad, and tll are linearly interpolated in time for simulation times between the times for which we have observed values. We then calculate the result of one minute of decay and diffusion on the updated expression values. Equations 1 (or 4) and 7 can be solved analytically if there is no production. We assume reflecting boundary conditions at the anterior boundary (the 35% line), because that boundary splits the anterior Hb and Gt peaks. Any protein diffusing out across the 35% line should be matched by protein diffusing in across the 35% line from the parts of the peaks not modeled. At the posterior boundary (the 92% line), we use an absorbing boundary condition (va(x,t) = 0), because that boundary does not intersect any of the gap gene peaks.

Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the National Science Foundation.