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

While it is encouraging that independent optimisation runs and methods give consistent parameter estimates, it is necessary to test the reliability and accuracy of these estimates using practical (or a posteriori) parameter identifiability analysis (see Introduction). We have performed such an analysis using two complementary approaches.

One approach to the practical analysis of parameter identifiability is based on a geometrical interpretation of the ‘optimisation landscape’ given by the value of the weighted-least-squares cost function (, in equation (2)) [35], [36]. To illustrate this approach, we will assume a two-dimensional parameter space for simplicity. This results in a three-dimensional topography of the optimisation landscape, where minima lie in ‘troughs’ or ‘depressions’ of the contour determined by the cost function. The more shallow the trough in which a minimum lies, the more uncertain the parameter estimate, since changing parameter values around the optimum will lead to only a slight increase in the value of the cost function. It is possible to characterise the local surface of any optimisation landscape around a given minimum using linear approximations. This allows us to define an ellipsoidal confidence region around our minimum, resulting in estimates for the confidence intervals for each of our parameters.

If there is no correlation among parameters, the principal axes of the confidence ellipsoid will lie parallel to those of parameter space. Confidence intervals for parameters can then be calculated as the intersect of the ellipsoid with these axes. Correlations among parameters are detectable as an inclination between the ellipsoid's principal axes and the axes of parameter space. This makes it possible to calculate two distinct ranges: the dependent confidence interval is given by the intersection of the ellipsoid with a given parameter axis (as above), while the projection of the ellipsoid region onto the parameter axis specifies the independent confidence interval. Independent confidence intervals typically overestimate the uncertainty in parameters, while dependent confidence intervals underestimate it. If both confidence intervals turn out to be similar and small, a parameter can be considered well determined.

Confidence intervals, as calculated by equations (4) and (5) (see Materials and Methods) are shown in Figure 4. Compared to the entire range of search space, confidence regions for all three rate parameters (, , and ; see equation (1)) are small in models of all three gap genes. Dependent and independent confidence intervals for and deviate significantly, suggesting strong mutual correlation among model parameters. This is not the case for , where independent and dependent confidence intervals are very similar. Note that the lower limits of some intervals for and are negative, and therefore lie outside the allowed range of parameter values. This artifact results from the linear approximation of the optimisation landscape used in this method. Confidence intervals for delay parameters are larger compared to rate parameters. Nevertheless, they lie within a well confined and biologically plausible range (Figure 4). As for the case of and above, there appears to be a high degree of correlation for with other parameters in all three models.

Correlation coefficients between parameters can be calculated from the covariance matrix (Figure 5A; see equation 6 in Materials and Methods). In all three models, correlation is high between and . This is expected since high decay rates can compensate for high production rates. Both of these parameters are also correlated to the delays given by . These correlations are highest for gt, and still very substantial for both Kr and kni. Again, this is to be expected since production delay can be mimicked to some degree by low production rates. In contrast, we found that diffusion rates are largely independent of other model parameters, except for a slight negative correlation between and , and between and . This could be due to the extremely low values of diffusion rates in all of our models, or due to the fact that diffusion affects spatial, rather than strictly local, regulatory mechanisms, which could explain the increased degree of decoupling between the two processes.

While computationally efficient, the linear identifiability analysis described above can lead to serious artifacts or biases in the estimation of confidence intervals due to its simplifying assumptions. Therefore, we validated its results by using the computationally much more expensive approach of bootstrapping [31], [71]–[73]. The bootstrap method is based on resampling protein expression patterns from distributions defined by the mean and variance of our measurements (for an equivalent analysis of the sensitivity of parameter estimates with regard to perturbations in the mRNA data see Supplementary Text S4). The model is then fitted to a large number ( in our case) of such sampled noisy patterns. Confidence intervals and correlations for parameter estimates can be directly extracted from the resulting parameter distributions.

Distributions of parameter estimates obtained by bootstrapping are shown in Figure 6. In all cases, estimated parameter values are confined to relatively small subregions of search space. Only diffusion rates show a tendency towards saturation at their lower limit (; Figure 6G–I). Distributions are generally unimodal, with the exception of Kr which shows two distinct clusters in parameter space. The cause of this bimodal distribution remains unclear. One of the clusters (568 solutions) has implausibly small values for ( min). Therefore, we only considered the remaining 432 solutions (see dashed circle in Figure 6D) for further analysis. Another interesting feature of these parameter distributions are the two horizontal lines visible in Figure 6E, which indicate an exclusion of values around 6.25 and 6.30 minutes. This corresponds roughly to the time between data points, which seems to indicate that the structure of the data used for model fitting has a non-negligible impact on parameter estimation in this case.

None of these irregularities observed in parameter distributions seriously affects our ability to compute confidence intervals. This was done by determining the 95-percentile range for each parameter separately. The resulting confidence intervals and the initial guess are shown in Figure 4. The optimal solution of the unperturbed data set in every case lies within or very close () to the limits of the corresponding confidence interval (diamonds in Figure 4). With the exception of confidence intervals around , the size of bootstrap intervals lies between those of the dependent and independent intervals calculated by linear approximation. In general, this confirms the accuracy and reliability of this method.

Some notable exceptions apply. First, most confidence intervals based on bootstrapping are clearly asymmetric around the estimated optimal values. This asymmetry reflects a non-linear dependence of the optimisation problem on parameter values, which cannot be captured by confidence intervals calculated from linear approximation. Second, size of confidence intervals for obtained by the bootstrap are frequently more similar to the size given by the dependent confidence interval, indicating that delay ranges are more accurately determinable than estimated by linear approximation.

Correlation matrices calculated by linear approximation or bootstrapping are very consistent. The anisotropic shape of parameter distributions resulting from bootstrapping reveal strong positive correlations between rate parameters for production and decay ( and ; Figure 5B, 6A–C). Somewhat weaker, but nevertheless strong correlations occur between these same rate parameters and the production delays (Figure 5B, 6D–F). In contrast, diffusion rates are much less correlated with any of the other parameters (Figure 5B, 6G–I).

In this study, we have used a reverse-engineering approach to test whether post-transcriptional regulation is required for the correct expression of gap protein domains. For this purpose, we have created a high-resolution quantitative dataset of mRNA expression patterns for the gap genes Kr, kni, and gt covering the entire blastoderm stage. Comparison of gap mRNA and protein expression data indicates that both are remarkably similar, although features in the mRNA data emerge a few minutes earlier than those of the corresponding protein patterns. Results of our model fits confirm this general impression: the timing and position of gap protein domains can be explained largely by a simple linear delay model, which assumes that protein patterns correspond to mRNA patterns a given amount of time in the past (plus a small contribution of protein diffusion). Based on this, we conclude that post-transcriptional regulation is not essential for gap gene mediated pattern formation. This result confirms a widely held assumption by the Drosophila research community that had never been put to a rigorous test.

On the other hand, our results reveal surprising and significant differences between mRNA and protein levels. In particular, our models fail to correctly reproduce both early dynamics of expression initiation (for Kni and Gt), and late maintenance of protein levels (for all three gap proteins; Figure 3; see also [58]). This indicates that temporally and spatially specific post-transcriptional regulation is required to explain these particular expression features. Only a detailed quantitative study, such as the one presented here, is able to detect such subtle nuances. No experimental evidence is currently available on the regulatory mechanisms or the functional importance of these newly discovered expression features.

Our models yield predictions that are informative and specific enough to enable focused molecular and biochemical investigations of these phenomena. The early boost in build up of Kni and Gt protein could be explained by a modulation of protein production rates. One potential mechanism for this would be a translation rate that depends on the diffusion limited arrival of mRNA molecules at the ribosomes. The maintenance of high protein levels despite rapid mRNA decay towards the end of the blastoderm stage indicates some mode of protein stability regulation. Such temporal regulation has been observed for Bcd protein [74], [75]. It could be achieved by non-linear dependence of the decay rate on protein concentration. Co-operative stability—increased longevity of dimers compared to monomers—has been proposed as a potential mechanism for increased protein half life at high concentration levels [76].

One last aspect of post-transcriptional regulation that requires our attention is the protein production delay predicted by our models. These delays, between and min long (Table 1), are short but yet significant enough to affect the dynamic regulatory properties of the system. They have several effects: First of all, production delays must be kept rather short to allow pattern formation on a time scale of less than min in a system of rapidly dividing nuclei [58], [59]. There is some experimental evidence to show that this is achieved through a compact gene structure—short open reading frames with a very limited number of short introns [77]. While kni has a primary transcript about 3 kilobases (kb) long, its paralogue knirps-related (knrl) (encoding a functionally equivalent protein) contains a long intron which results in a primary RNA of about 23 kb. Its limited length allows kni to become expressed early, at cleavage cycle 13. In contrast, cytoplasmic mRNA of knrl only appears around mid cleavage cycle 14A, about half an hour later. The second aspect of the production delay is important in the context of the transient nature of gap gene patterning. While it has been shown that mRNA and protein levels of a gene converge at steady state [78], they can be significantly different when a system is far from asymptotic behavior. In the case of the gap gene system, this is reflected by the systematic anterior displacement of mRNA compared to protein domains, caused by the dynamic anterior shift of gap gene expression domains over time [22], [25], [59]. This phenomenon had been attributed to post-transcriptional regulation by some authors [67], but can now be fully explained by a combination of the anterior movement of the domains, the production delay, and the slightly different half lives of mRNA and protein.

Finally, production delays that are on the same order of magnitude as the time scale of pattern formation can lead to severe alterations of the transient dynamic behavior of the system. For example, delays can greatly increase the time it takes for the system to reach its steady stage [79], [80]. This may be functionally important for gap gene patterning, where the expression domains in the posterior half of the embryo have to be kept moving anteriorly until the onset of gastrulation (and the subsequent disappearance of gap expression), while gap domains in the central part of the embryo remain stable and reach their steady states much earlier [25], [59].

As in the case of delays, our models yield predictions of rate parameter values that are plausible, informative, and experimentally testable. Predicted decay rates imply gap protein half lives that lie between and min, which is somewhat lower than the to min measured for the Bcd protein [74], [75]. Our predictions of diffusion and production rates are harder to assess. The reason for this is that they are formulated in relative units, since our measurements of mRNA and protein concentrations are relative and do not yield absolute concentrations. This limitation could be overcome by emerging experimental techniques that allow the estimation of absolute levels of mRNA and protein in vivo [61], [81]–[83]. However, we can already draw some conclusions from our estimated relative values. In particular, our results indicate that gap protein diffusion must be severely restricted. This is corroborated by our observation that gap protein domains are generally only about 1–2 nuclei wider than their corresponding mRNA domains (see Supplementary Table S1), and is consistent with the model-based prediction that diffusion is not required for correct gap protein mediated patterning [22], [25]. Finally, production rates are the most difficult to measure. In this regard, the prediction of our models that protein production and decay must be quite tightly balanced may be helpful to overcome this technical limitation.

Our simple model of transcriptional regulation is limited in several important ways. We have explicitly refrained from implementing particular post-transcriptional regulatory mechanisms due to the absence of specific experimental evidence at this point. Our main aim in this current study was to first establish whether any post-transcriptional regulation is necessary for gap gene regulation. Our results clearly show that such regulation is required for the proper level, but not timing and position, of gap gene expression. Future investigations will combine experimental and data-driven modeling approaches to extend the model, and render it more mechanistically accurate.

Another limitation concerns the coarse-grained nature of our production delay. It summarizes contributions by transcriptional elongation, mRNA processing and splicing, nuclear export, and translation (see Introduction). Such coarse-graining is warranted in the light of our experimental protocols, which cannot yet distinguish between nascent and primary transcripts, or mature mRNA. A more accurate measurement of mRNA production could be achieved by using intronic probes [84]. However, this is challenging in our context due to the extremely compact structure of gap genes. Another way to address this issue would be to formulate a model with a distributed production delay, reflecting the fact that the measured mRNAs are at different stages of their maturation process. We have not implemented such distributed effects in our current model since the benefit in terms of biological insight would be limited, while estimating a distribution of parameter values would pose significant technical challenges for model fitting.

The last, and most important, limitation of our current approach is that transcriptional and post-transcriptional regulatory processes involved in gap gene patterning are still implemented in different models. It is our aim to synthesise both of these stages into a regulatory network model featuring realistic production delays. We expect that such a model would solve several important issues. For instance, current gene network models still fail to reproduce the early regulatory and expression dynamics—based on regulatory inputs from maternal gradients only—in an accurate and biologically plausible manner [22], [36], [58]. An integrated model featuring a true expression delay will undoubtedly help to overcome this problem. Both the quantitative dataset of mRNA expression, and the mathematical analysis of gap gene translation presented here will be crucial for the implementation of such a model.

At a more general scale, we have provided a proof of principle that rigorous model fitting and parameter identifiability analysis are possible in the context of the complex regulation of animal development. We hope that this will enable a more widespread and rigorous application of reverse-engineering approaches to problems of biological pattern formation. In our view, this constitutes an important methodological advance, which is crucial to apply the considerable potential of quantitative reverse-engineering strategies for our understanding of development.

Blastoderm stage embryos of Drosophila melanogaster (raised at 25°C) were collected 1–4 hrs after egg laying. Embryos were fixed and stained using FITC- (Kr, gt) or DIG-labeled (kni) riboprobes, plus polyclonal antiserum against Even-Skipped (Eve) [85], according to standard experimental protocols [58], [86]–[88]. Nuclei were counterstained using Hoechst 34580. Imaging took place on a Leica TCS SP5 confocal microscope using a 20× objective, and an additional digital zoom of 1.3x. We imaged the blastodermal nuclear layer of laterally oriented embryos at two -positions, 1.0–1.2 μm apart. Data channels were scanned sequentially at a resolution of 1024×1024 pixels. Only embryos at cleavage cycle 14A (C14A) [89] were chosen for further processing. For earlier time points, we use previously published gap mRNA expression data [58].

Data processing and quantification methods are described elsewhere in detail [58], [69]. In brief, we create a binary whole-embryo mask by thresholding, which is used to automatically crop and align embryo images such that anterior is left, dorsal up. We identify nuclei and their surrounding territories of cytoplasm using watershed-based image segmentation algorithms [58], [57]. From these watershed masks, we determine the position of nuclei, as well as the concentration of mRNA (in nuclei plus surrounding cytoplasm), and Eve protein (in nuclei only). Embryos are classified into eight time classes (T1–8; each min long) during C14A, based on Eve expression patterns and morphological markers [69]. Non-specific background staining is removed as previously described [90]. Expression data are registered using a spline-based approach [91]. Background removal and data registration are implemented in an integrated tool [92]. Finally, registered mRNA data within a lateral strip (covering 10% of the embryo) are placed into 100 bins along the A–P axis, and concentration values in each bin are averaged per gene and time class. Individual expression profiles, integrated patterns, and the number of embryos used for each time class, are shown in Supplementary Figures S1, S2, S3. Gap mRNA expression patterns for C10–C13 were taken from [58], and scaled to provide a smooth transition between the two datasets. Integrated protein expression data used in our analysis are from the FlyEx database: http://urchin.spbcas.ru/flyex [57], [60]. We normalise our mRNA data (using the same scaling factor for all time classes) by adjusting peak concentrations to the maximum expression level observed for protein. Expression peaks, domain boundary positions (points of 50% maximum expression), and domain widths were calculated from spline approximations to the expression data as described in [63].

The basic objects of our model represent nuclei plus their associated surrounding cytoplasm (energids). The state variables of the model describe the concentration of intra-nuclear gap protein within each energid. Change in gap protein concentration across time and space is described by a system of ordinary differential equations (ODEs; see equation 1 in the Introduction), and depends on protein production from mRNA (concentration averaged across both nuclear and cytoplasmic portions of the energid), protein diffusion between energids, and protein decay. The model spans the entire blastoderm stage, from min after the onset of cleavage cycle 10 (C10; min) to the end of C14A at the onset of gastrulation ( min) [89]. During this time, four mitotic divisions occur (division 10, min; division 11, min; division 12, min; division 13, min; [89]). During mitosis, transcription of mRNA is interrupted and unfinished transcripts are actively degraded [93]. This process has never been quantified. Here, we assume fast mitotic mRNA degradation: therefore, the protein production term in our model (see equation 1) is set to zero, whenever the time point (current time minus the production delay) comes to lie within the time of mitosis. At the end of each mitotic phase, nuclei divide instantaneously (and thus the number of ODEs in the model increases approximately two-fold), and the distance between them is halved. Due to the presence of diffusion, our ODEs are coupled across space. The spatial range of our model is defined for each gap gene independently. In general, models cover most of the segmented trunk region of the embryo. Ranges were defined to include the posterior boundary of the anterior Gt domain, the central Kr domain, the abdominal Kni domain, and the posterior domain of Gt (Kr: 25.5–88.5%, kni: 32.5–88.5%, gt: 32.5–95.5% A–P position, where 0% is the anterior pole).

Although our models are feed-forward and linear, they depend on a non-linear external input (the mRNA expression profiles). Therefore, we solve these systems of ODEs numerically using an implementation of the MATLAB dde23 solver in C [94]. The solver was modified to satisfy the following requirements: (i) it must be able to provide dense output at time points for which there is no external input data; (ii) it must be able to handle discontinuities propagating through the system due to the delay history; and (iii) it has to handle implicit formulas if the stepsize becomes bigger than the delay. We use linear interpolation between data points to provide mRNA concentrations as external inputs at arbitrary points.

Structural parameter identifiability analysis was performed using the Laplace transform based approach [28]. The idea is to obtain the transfer matrix of the system in rational canonical form and to assess whether the transfer matrix is unique. Supplementary Text S2 provides a detailed description of this approach, and the calculations we performed.

Numerical solutions are produced for time points C10–C13, and T1–T8 within cleavage cycle 14A. We then calculate a weighted sum of squared differences (according to equation 2), which is minimised using two alternative optimisation strategies. The first of these consists of global optimisation using (parallel) Lam Simulated Annealing (pLSA; [95]–[98]). This method is reliable and robust, and has been successfully used in previous reverse-engineering studies of the gap gene system [22], [58], [64]. pLSA is computationally intensive and was implemented in C. The second approach consists of a scatter search approach (eSS), which systematically explores parameter space, triggering a local gradient-based search [99] whenever a promising parameter set has been found [100]–[104]. eSS is implemented in the AMIGO toolbox (based on MATLAB; [105]). Both strategies resulted in virtually identical model fits and parameter estimates (Figure 3; Table 1).

The following search space limits were used for optimisation. : 0.005–5.0; : 0.0347–0.6931; : 0.0–0.3; : 0.0–10.0.

The quality of the fit between data and model output is measured by the root mean square () score, which represents the average difference between modeled and measured protein concentrations across all data points:(3)where denotes the concentrations of protein in nucleus at timepoint , while corresponds to the according simulated value of protein concentration dependent on the chosen parameter set . The RMS—unlike the weighted least square sum (see equation (2))—is independent of not only the noise in measurement but also the number of data points used for fitting. It therefore makes model fits of different gap genes comparable to each other on a quantitative basis.

We used two alternative strategies for practical parameter identifiability analysis. The first one is based on a local linear approximation of the ‘energy’-landscape given by the objective function (2) around a given optimum as described in [35], [36], [66]. Estimates of the local contour of this landscape are used to determine an ellipsoidal confidence region around the optimal parameter set obtained by eSS. Dependent confidence intervals () are then given by the intersection of the ellipsoid with the parameter axes(4)while independent confidence intervals are specified by the projection of the ellipsoid onto the parameter axes(5)

Correlations among model parameters can be calculated based on the covariance matrix as(6)In all equations, diagonal entries in correspond to with beeing the standard deviation on measurement. is the inferred parameter set obtained via eSS and denotes the sensitivity matrix of the model given by the first order derivative of the observables with respect to the parameters. is calculated as the upper part of Fishers distribution with and degrees of freedom(7)with(8)