Reverse-engineering post-transcriptional regulation of gap genes in Drosophila melanogaster

Systems biology proceeds through repeated cycles of experiment and modeling. One way to implement this is reverse engineering, where models are fit to data to infer and analyse regulatory mechanisms. This requires rigorous methods to determine whether model parameters can be properly identified. Applying such methods in a complex biological context remains challenging. We use reverse engineering to study post-transcriptional regulation in pattern formation. As a case study, we analyse expression of the gap genes Krüppel, knirps, and giant in Drosophila melanogaster. We use detailed, quantitative datasets of gap gene mRNA and protein expression to solve and fit a model of post-transcriptional regulation, and establish its structural and practical identifiability. Our results demonstrate that post-transcriptional regulation is not required for patterning in this system, but is necessary for proper control of protein levels. Our work demonstrates that the uniqueness and specificity of a fitted model can be rigorously determined in the context of spatio-temporal pattern formation. This greatly increases the potential of reverse engineering for the study of development and other, similarly complex, biological processes.

Systems biology is characterised by the tight integration of experiments and computational modeling. One way to achieve such integration is through reverse-engineering approaches, where dynamical models of regulatory or biochemical reaction networks are fit to quantitative data [1]–[9]. Reverse engineering has been successfully used for systems analysis in many contexts, from microbial metabolic, signaling and regulatory networks (see, for example, [10]–[20]) to pattern-forming developmental processes in animals (e.g. [21]–[25]). The approach is illustrated by the systems biology modeling cycle shown in Figure 1 (adapted from [26]). As a first step, a mathematical model is formulated that incorporates the basic assumptions and hypotheses we have about the regulatory process under study. The model is then tested by fitting it to metabolic or expression data. This is achieved by repeatedly altering its parameters and selecting suitable (mainly better-fitting) solutions. A successful fit will yield a unique set of parameter estimates that cause the model to reproduce the data accurately. In this case, model output and estimated parameter values can be analysed to gain biological insight. For instance, regulatory parameters contain information on the strength and type of regulatory interactions in a network. If the model fails to produce a unique solution—predicting a large set of variant networks instead—it is underdetermined and more data need to be collected. If the model cannot fit the data, the underlying hypothesis needs to be adjusted, or additional mechanisms and factors need to be incorporated. Successive model-fitting/data-acquisition cycles yield an increasingly accurate quantitative picture of the underlying regulatory network.

While this approach has great potential for the investigation of complex biological regulatory systems (e.g. [27]), it also harbors many significant and non-trivial challenges. One of those is that it is often difficult to decide what kind of data, and what kind of model are needed to enable a successful fit. Another challenge is to analyse whether a given solution is indeed specific and reliable. There are a number of mathematical methods designed to establish whether a reverse-engineering problem is well posed—in other words, whether it is able to produce a unique and consistent solution [5], [7], [28]–[30]. First, structural (or a priori) parameter identifiability analysis can be used to examine whether the problem has a non-trivial solution at all [31]–[33]. Second, practical (or a posteriori) parameter identifiability analysis tells us whether estimated parameter values are significant and reliable [31], [34]–[36]. Finally, methods for optimal experimental design are employed to determine what kind of measurements (for which regulatory factors and which time points, for example) would improve the quality of the fit most significantly [5], [37]–[40].

So far, unfortunately, the application of these powerful methods to gain specific and novel biological insights has been limited. This is due both to the complexity of most real-world biological regulatory systems and the nature of the data used in reverse-engineering studies. Most of these studies use models based on large systems of coupled non-linear differential equations. This makes it challenging to apply structural identifiability analysis. Moreover, model fitting is generally computationally intensive due to the significant number of parameters to be estimated. This renders rigorous practical identifiability analysis extremely time consuming. And finally, high-throughput datasets used for model fitting often exhibit high levels of measurement error, combined with low numbers of replicates. Under these circumstances, it is difficult to accurately assess data variance, which is required for both practical identifiability analysis and optimal experimental design. For all these reasons, reverse-engineering studies often proceed on an empirical basis, without being able to rigorously establish parameter identifiability or the suitability of the datasets and models used.

Here, we present a reverse-engineering study which combines model fitting by global optimisation strategies with rigorous structural and practical identifiability analysis. We apply this methodology to a complex regulatory problem: the dynamics of spatio-temporal pattern formation in the early embryo of the vinegar fly Drosophila melanogaster. The biological question we are addressing is the importance of post-transcriptional regulation in animal development. While many studies of pattern formation focus on differential transcriptional regulation of genes (e.g. [41], [42]), other levels of expression control—such as regulated RNA splicing, processing, translational regulation, or regulated stability and degradation of gene products—cannot be ignored [43]. There is increasing evidence that protein levels do not generally match those of their respective mRNAs [44]–[46], and many protein expression patterns do not even coincide with the timing and localisation of mRNA transcription [47], [48]. These discrepancies are due (at least in part) to control at the level of protein translation. Indeed, some of the earliest studies of translational control were carried out in Drosophila (reviewed in [48]). A number of pioneering studies examined the effect of translational repression on maternal morphogen gradients, such as those formed by the protein products of the maternal genes hunchback (hb) and caudal (cad). mRNAs derived from those genes are distributed uniformly while their proteins form steep concentration gradients with antero-posterior polarity [49]–[54]. More recently, systems-level studies indicate that such post-transcriptional regulation is widespread and of general importance.

Protein expression levels in yeast cannot be predicted from mRNA concentrations alone [55], and a similar lack of correlation between mRNA and protein is observed in systems as different as the minimal bacterium Mycoplasma pneumoniae [47] and mammalian cell lines [46]. Therefore, post-transcriptional regulatory mechanisms must be incorporated in a systems-level understanding of cellular and organismal function.

In this study, we investigate the role of post-transcriptional regulation within the context of a well established experimental model system: the gap genes involved in segment determination during the blastoderm stage of early Drosophila development (reviewed in [56]). Since the relevance of post-transcriptional regulation for maternal hb expression is well established (see above, and [52], [53]), we will focus on the remaining three trunk gap genes: Krüppel (Kr), knirps (kni), and giant (gt). All these genes encode transcription factors, and are expressed in broad, overlapping domains along the embryo's antero-posterior (A–P) axis. Gap genes respond to activating transcriptional regulatory inputs from long-range maternal morphogen gradients—such as Bicoid (Bcd), Hb, and Caudal (Cad)—as well as repressive inputs from the terminal gap genes tailless (tll) and huckebein (hkb). In addition, there is extensive repressive gap gene cross-regulation, which is required for the correct dynamic positioning, maintenance, and sharpening of each gap gene expression domain.

The advantages of using the gap gene network for our case study are twofold. The first advantage is that gap gene patterning is relatively simple and tractable compared to other developmental processes. It essentially occurs in one dimension, along the A–P axis of the embryo. No significant tissue rearrangements or growth are involved. Diffusion is not yet limited by cell membranes as the embryo is still syncytial at this stage. In addition, all three genes considered here have a very compact structure, with only one or two short introns, and none of them exhibits any sign of alternative splicing. The second advantage is that the gap gene system is exceptionally well understood. All genes involved in segment determination have been identified and their interactions have been characterised at the genetic and molecular level (see [56] and references therein). More importantly, there exist extensive quantitative datasets (including accurate variance measurements) for spatio-temporal gap mRNA and protein expression [56]–[63]. These datasets have been used to fit a range of gene regulatory network models, analysis of which has led to many quantitative systems-level insights into the dynamic mechanisms underlying gap gene regulation [22]–[25], [36], [58], [64]–[66].

All these previously published models focus on transcriptional regulation of gap genes. They lump together transcriptional and post-transcriptional phases of gene regulation, and take into account only protein concentrations (not mRNA) as model observables. Therefore, these models implicitly assume that post-transcriptional regulation is not required to explain the patterns formed by gap genes. This assumption is not unreasonable, given the similarity of gap mRNA and protein patterns, and the fact that such simplified models can reproduce gap protein patterns to a high degree of accuracy and temporal resolution [22], [36]. Moreover, the experimental literature contains very little evidence or arguments for post-transcriptional regulation of Kr, kni, or gt. The only exception we could find is a paper by Gaul et al. [67], which invoked post-transcriptional regulation of Kr to explain the anterior displacement of its mRNA domain with respect to the protein pattern. This phenomenon was later shown to be due to the dynamic anterior shift of the central Kr domain [22]. However, absence of evidence is not evidence of absence. Therefore, it is necessary to put the hypothesis that post-transcriptional regulation is not required for gap gene patterning to a rigorous and quantitative test.

As mentioned above, we test this hypothesis using a reverse-engineering approach. This is achieved by formulating a model, which incorporates the simple assumption that gap protein patterns reflect those of their respective mRNAs a given amount of time earlier in development (plus a small contribution by gap protein diffusion; see Figure 1, right-hand panel). Since we do not consider gap gene cross-regulation, we can model Kr, kni, and gt separately. Each model is then fitted to protein expression data, using mRNA patterns as external inputs, or boundary conditions. If our models are able to reproduce gap protein patterns correctly, we can conclude that no post-transcriptional regulation is required for the expression of the gap genes considered here. If our models fail to fit, however, we will be able to identify those expression features that do rely on post-transcriptional regulatory processes.

Our model consists of the following system of ordinary differential equations, representing the change of gap protein concentration over time in a row of dividing nuclei along the A–P axis of the embryo:(1)Here, the dependence of protein on mRNA patterns is linear. is a time-delayed external input representing mRNA concentration minutes ago. We obtain for arbitrary time points by linear interpolation from measured expression data points. Transcription is paused during mitosis, when chromosomes are condensed. Therefore, the production term is set to zero whenever falls into a mitotic period (see Materials and Methods for time schedule). is a production delay which summarises the time from initiation of transcription (when transcripts become detectable by our staining methods) to the appearance of the resulting protein. It includes contributions from the completion of transcription, RNA splicing and processing, nuclear export, and translation. is the rate of protein production from mRNA. is the rate of gap protein diffusion among neighboring nuclei , , and . It depends on the square of the distance between nuclei which is halved upon each nuclear division [64]. is the protein degradation rate.

After establishing the structural identifiability of our parameters, we need to estimate (production delay), (production rate), (diffusion rate), and (degradation rate) by fitting the model to data. This is achieved by minimising the weighted sum of squared differences between experimentally measured () and modeled protein expression values () for each nucleus and time point for which we have data:(2) is the parameter vector to be estimated. Weights are given by measured variances for each data point. Equation (2) represents a weighted least squares (WLS) problem, which we solve using global optimisation methods as described below. The reliability and accuracy of the resulting parameter estimates are then analysed using practical parameter identifiability analysis.

Using the reverse-engineering approach described above, we have obtained fitted models and parameter estimates for post-transcriptional regulation of Kr, kni, and gt. Identifiability analysis shows that our fitting results and parameter estimates are robust and specific. They yield values for rate parameters which are biologically plausible and informative with regard to the time scale and diffusive properties of gap gene patterning. Our fits reveal that post-transcriptional regulation is not required for the correct timing and positioning of gap protein domain boundaries. They do suggest, however, that post-transcriptional regulation is required for the accurate control of gap protein levels, implying some temporal regulation of translational efficiency, and/or protein stability. Specifically, our models predict an early boost in translational efficiency, plus a general stabilisation of gap protein products towards the end of the blastoderm stage.

We have created a quantitative mRNA expression dataset with high spatial and temporal resolution for the trunk gap genes Krüppel (Kr), knirps (kni) and giant (gt), which spans the entire duration of the blastoderm stage in the early embryo of Drosophila melanogaster (cleavage cycles C10–C14A; C14A is further subdivided into time classes T1–8). In contrast to previously published semi-quantitative gap gene mRNA data—based on colorimetric (enzymatic) staining protocols, wide-field microscopy, and an efficient but simple data processing pipeline [66], [68] —we used fluorescent staining protocols, confocal microscopy and fully quantitative data processing methods (see Materials and Methods, and [69]). This work extends a previously published fully quantitative expression dataset for gap gene mRNAs, which only covered the early part of the blastoderm stage (C10–C13) [58]. Our data consist of quantified time-series of gap gene mRNA expression (Figure 2; Supplementary Figures S1, S2, S3), which are equivalent and comparable in quality, as well as spatio-temporal range and resolution, to the comprehensive protein expression data available from the FlyEx database (http://urchin.spbcas.ru/flyex, [57], [59], [60]). This allows us, for the first time, to rigorously and accurately compare gap gene expression during the blastoderm stage at both mRNA and protein level.

Such a comparison between spatio-temporal gap mRNA and protein expression reveals that, to a first approximation, the mRNA patterns appear very similar to those observed for the corresponding proteins: all mRNA transcripts are expressed in broad, overlapping domains, whose relative timing and spatial arrangement with regard to each other strongly resemble that of gap protein domains (Figure 2A). In addition, the central, abdominal, and posterior mRNA domains of Kr, kni, and gt shift towards the anterior of the embryo over time (Figure 2B). The extent of this movement is on the same order of magnitude as the analogous shifts of the corresponding protein domains (Figure 2B, Supplementary Table S1) [22], [59], [66]. Since anterior domain boundaries generally move less far than posterior ones, domain width of both mRNA and protein domains decreases over time (Figure 2B, Supplementary Table S1) [59]. Based on these observations, we can conclude that expression dynamics of gap mRNA and protein domains are largely qualitatively equivalent with regard to each other.

However, if we examine the data more closely, two significant differences between mRNA and protein become apparent. First, boundary positions of mRNA domains—with the exception of the anterior Kr border—are displaced anteriorly compared to their corresponding protein domains (Figure 2B; Supplementary Table S1). This displacement is caused by the anterior shift of gap domains over time [22], [59]. The effect is substantially more pronounced for posterior domain borders than for anterior ones. It is associated with a generally slightly smaller width of mRNA domains compared to protein (Supplementary Table S1). Second, the timing of initial and maximum peak expression is delayed for protein compared to mRNA (Figure 2C). Delayed first appearance of protein versus mRNA patterns during the early blastoderm stage has been reported and quantified previously [58], [59]. Our data reveal a similar phenomenon in late-blastoderm expression dynamics: mRNA expression of all three gap genes peaks around 30 min before gastrulation (time class T3 in Figure 2C), while protein expression shows a maximum approximately 15–20 min later (time class T6–7). This obviously agrees with the fact that it takes time to export the mRNA from the nucleus, to process and translate it into protein. In addition, we detect a post-peak decrease in mRNA abundance that was not reported in an earlier PCR- and microarray-based analysis of pre-gastrulation gene expression with a lower temporal resolution than our data [70]. Interestingly, this trend is not reflected in levels of gap proteins, which only show a marginal decrease (if any) before the onset of gastrulation (see Figure 2C and also Supplementary Text S1).

In addition to our quantification of the timing and position of averaged gene expression patterns, we have also analysed the embryo-to-embryo variability of gap domain width as well as peak and boundary positions during the blastoderm stage. This had not been possible with our earlier, semi-quantitative dataset [66]. It has been reported previously that the precision of gap protein domain boundary positions increases over time due to cross-regulatory interactions [24], [25], [59]. Such a trend—although much less obvious—is also present in our mRNA data (Supplementary Figure S4). Noise levels in the mRNA data are generally higher than in protein patterns. Both Kr and kni (but not gt) show higher levels of variability in mRNA compared to protein data at a majority of sampled data points. Moreover, we observe a high level of fluctuations between time points in our mRNA data (for all three genes) indicating increased levels of experimental noise. This may be due to the harsher treatment of embryos for in situ hybridisation compared to antibody staining (see [58] for a more detailed discussion). Nevertheless, the overall trends are similar for mRNA and protein, as levels in variability of all measured expression features is generally lower at T6 than at earlier time points during C14 (Supplementary Figure S4). This indicates canalisation of development at two levels: first, protein patterns are generally more precise than mRNA domains at comparable stages, and second, mRNA precision—as is the case for protein—increases over time.

The fact that mRNA and protein patterns of the gap genes Kr, kni, and gt are similar (yet not identical), and show a delay in dynamics with regard to each other, raises the non-trivial question whether protein patterns simply reflect earlier mRNA levels (with a small additional contribution by protein diffusion), or whether spatially and temporally specific post-transcriptional regulation is required to account for the observed distribution of proteins. We use a reverse-engineering approach to distinguish between these two alternative possibilities. In this approach, we test the hypothesis that no post-transcriptional regulation is required by fitting a simple dynamical model to data. This model incorporates the following assumptions (see Introduction): It includes time-delayed but linear production of protein from its mRNA, as well as protein diffusion and decay. It takes mRNA patterns as an input to predict protein expression. A good fit of the model to protein expression data would therefore favor absence of post-transcriptional regulation, while a failure to fit the data would point us to specific features of gap protein expression that require regulated nuclear export, splicing, or translation.

Our reverse-engineering approach can only give us quantitative and specific insights into the problem of post-transcriptional regulation if it is fit to data in a manner which is as rigorous and reproducible as possible. As a first step, this requires us to determine whether the model is formulated in a way such that the fitting procedure has a unique solution. Since our model is feed-forward and linear, this can be achieved using structural (or a priori) parameter identifiability analysis. This analysis is performed under an ideal scenario where noise-free time-continuous experimental data are assumed to be available, and the objective is to answer the question whether under those ideal conditions the parameters can be given unique values. There are three possible outcomes: (1) The model is structurally globally identifiable (s.g.i.) if all parameters can be uniquely identified within a biologically meaningful region of parameter space, which we will call (). (2) The model is structurally locally identifiable (s.l.i.) if one or more parameters can be uniquely identified in a given neighborhood , or (3) the model is not structurally identifiable, if neither of (1) or (2) apply.

Although several methods for the analysis of structural identifiability of linear models exist, the model described by equation (1) presents the peculiar challenge of incorporating a delay parameter within the input function (the production term). In this scenario the Laplace transform () based method may be used to assess whether the s.g.i. condition holds (see [28] and references cited therein). The underlying idea is to verify whether a canonical form of the transfer matrix of the system 1 is unique. The basic steps are the following: the model (equation 1) is rewriten in matrix form and its Laplace transform is computed; the possibility of computing the transfer matrix is demonstrated by an invertibility conditon; the analytical canonical form of the transfer matrix is then obtained; symbolic manipulation is finally used to prove uniqueness of the transfer matrix. Details of our calculations can be found in Supplementary Text S2.

In the case of our model for post-transcriptional regulation (see Introduction, equation 1) structural identifiability analysis reveals that model parameters are globally identifiable in the realm of non-negative real numbers (which is expected given that rate and delay parameters cannot be less than zero). For this result to apply, the following conditions must be met by the experimental data: (1) some concentrations of mRNA and proteins must be non-zero in the interior of the model range at sampling times, and (2) protein data must be available for time points before or after mitosis (i.e. when the production rate is not zero; see Supplementary Text S2 for details). Both of these conditions are met in our data, and we conclude that our model parameters are globally structurally identifiable, i. e. our optimisation problem has a specific and unique solution.

Next, we proceeded to fit our model to quantitative spatio-temporal protein expression data. This was done for Kr, kni, and gt separately by solving the model for each gene numerically, and minimising a weighted sum of squared differences between expression patterns predicted by the model and those measured by experiment (see equation 2). Minimisation of squared differences was achieved by two different global optimisation algorithms, based on parallel simulated annealing (pLSA), and an enhanced scatter search (eSS) method, respectively (see Materials and Methods for details and references). Both of these independent optimisation approaches resulted in equivalent model fits and parameter estimates (Figure 3; Table 1). In order to further corroborate the robust performance of our algorithms, we performed a systematic comparison of eSS with a number of standard global optimisation approaches. The results can be found in Supplementary Text S3. They indicate that different algorithms show significant differences in computational efficiency, but converge to very similar solutions.

Fitting results differ slightly between genes. The best fit between model and data is obtained by Kr optimisation runs, with a root-mean-square (RMS) score of around 11.1–11.2 (Table 1; see Methods for a mathematical definition of RMS, which represents the average deviation of model from data for each data point). Although minor patterning defects can be observed at early stages (especially between C13 and T2) and expression levels disagree somewhat at the last time points (T7/T8, see below), model and data match to within the noise level of the data at intermediate times (Figure 3, left column).

Of all fitting solutions, Kni shows the largest overall deviation between model and data with a RMS score of 21.1–21.2 (Table 1). Nevertheless, position and shape of the Kni protein domain, as well as the temporal dynamics of expression, are reproduced correctly (Figure 3, centre column). In particular, the model shows peak expression at the correct stage, and agrees with data to high accuracy in non-expressing areas. In contrast, protein expression predicted by the model is consistently lower than measured within the area of the abdominal Kni domain. This discrepancy accounts for the large RMS value.

Expression patterns resulting from Gt optimization runs exhibit similar properties as those of Kni, with a slightly lower residual error (RMS around 17.2; see Table 1). Again, the timing, position, and shape of both Gt domains are reproduced quite accurately in the model, while predicted expression levels are generally too low (Figure 3, right column). The only noticeable positional defect is a slight anterior displacement of the posterior gt domain at early stages (up to T1).

Despite the problems with reproducing accurate levels of expression for Kni and Gt, all three models show initiation and build-up of gap proteins at the appropriate stages of development, and the qualitative shape of the temporal expression profile is reproduced correctly up to T6 (Figure 3B). In contrast, model and data disagree conspicuously at later time points (T7/8), since the model predicts rapidly diminishing concentrations of gap proteins before the onset of gastrulation. This downregulation is much weaker (Kni), or entirely absent (Kr, Gt) in the data.

In summary, our models capture the position, shape, and width of gap protein domains accurately. Minor deviations in these spatial expression features are only observed during earlier time points, when noise levels in the data are high. Temporal features of gap protein expression—such as initiation of expression, shifts in domain position, or the time point of maximum expression—are also reproduced correctly. However, our models fail to reproduce the exact levels of gap protein expression, as well as their downregulation towards the end of C14. These two specific failures of our model fits have interesting implications for our understanding of gap gene regulation (see Discussion).

Model fitting resulted in reproducible and biologically plausible estimates of parameter values for Kr, kni, and gt. We performed a number of independent optimisation runs for each model using both of our two alternative fitting strategies (100 runs of pLSA, 10 runs of eSS per model). Parameter estimates from different runs varied only minimally between solutions (for pLSA see Supplementary Figure S5), and estimated parameter values provided by either of the two alternative optimisation strategies agreed to high accuracy. This indicates that our parameter estimates are robust with regard to the choice of optimisation strategy. Table 1 shows parameter values from a representative optimisation run for each fitting approach.

Predicted values for the delay parameter (see equation 1) require some more detailed attention, since such parameters are notoriously difficult to estimate. For this reason, we verified the validity of our estimates for by the following numerical approach: we performed a series of fits for the gt model, with fixed to values between 0 and 8 minutes, including a particularly high-density sampling of around the values predicted by optimisation (Supplementary Figure S6). These control fits show a minimum of the residual error which coincides precisely with the parameter values inferred from optimisation. This corroborates the reproducibility and accuracy of our approach.

Our parameter estimates are informative from a biological point of view, and yield experimentally testable predictions. First of all, we note that production and decay rates ( and , respectively) are well balanced in this system. Decay rates correspond to protein half lives of 9.1 (Kr), 9.0 (Kni), and 6.2 (Gt) minutes, indicating that gap proteins must be very unstable—to enable patterning at the extremely short time scale of gap gene expression dynamics. Gap protein diffusion is generally very low, especially in the case of Kr. Protein production delays (incorporating contributions of transcription, splicing, nuclear export, and translation; see Introduction) range between and minutes. While the upper value is within the expected range (see Discussion) the former estimate is rather low, and may need further investigation (see also next section).