Efficient reverse-engineering of a developmental gene regulatory network
Statistical analysis of parameter estimates was performed as described in [49], [64]. Here we only provide a short description of the calculation of parameter confidence intervals.
Reverse engineering results in a vector of estimates . Once the parameter vector is found, a posteriori identifiability analysis [65]–[67] reveals how reliable the obtained estimate is. The ellipsoidal region around in which the true parameter vector lies with a certain probability (we set ) is defined by(6)where(7)and and the number of parameters and data points, respectively. is the upper part of Fisher's distribution with and degrees of freedom. J is the Jacobian matrix of size defined as , where is the vector of weighted discrepancies between model output and data. From equation 6 one can derive dependent and independent confidence intervals for each parameter estimate . These are, respectively,(8)and(9)Here and are obtained from the singular value decomposition of . It is well known that in the presence of strong correlations between parameters, the dependent confidence intervals underestimate the confidence region while the independent confidence intervals overestimate it. For detailed explanations of these statistical quantities and their derivations we refer the reader to [49], [64], and references therein.
Image processing and extraction of expression domain boundaries were performed using a custom-made processing pipeline with a graphical user interface developed in Java (using the ImageJ API; http://rsbweb.nih.gov/ij). Intermediate processing steps and domain boundary positions were stored in a MySQL database, with a web interface (SuperFly) developed by the CRG Bioinformatics Core Facility. SuperFly is available online at: http://superfly.crg.es. We used scripts written in Python, Perl and R for the preparation of integrated data sets, for the generation of artificial external inputs, and for analysis of gene circuit models. Code for numerical solution and optimisation of gene circuits by pLSA [39], [45]–[49], [51], [62] is implemented in C, using the GNU Scientific Library (GSL, www.gnu.org/software/gsl), the Sundials ODE Solver Library [68], and the Open MPI message-passing interface (www.open-mpi.org). For numerical integration of ODEs, we use an implicit variable-order, adaptive-stepsize, multi-step method; a band-direct solver calculates the set of equations that is generated at each integration step [68].
Over the last two decades, the potential of reverse engineering has been demonstrated by a pioneering case study—led by John Reinitz and colleagues—where gene circuits have been used to characterise and analyse the gap gene network in Drosophila melanogaster [40], [45]–[51]. Despite this, the gene circuit method has not yet been applied more widely. One reason for this is that it took many years to establish the required quantitative data set of spatial gap protein expression patterns [41]–[43].
We have developed a simplified protocol for data acquisition and processing, which allows us to create a quantitative data set of spatial gene expression patterns in a time span of months rather than years. Instead of using protein expression data, we have quantified mRNA expression patterns by colorimetric (enzymatic) in situ hybridisation, imaged using a wide-field, compound fluorescence microscope. The resulting data set is of reduced quality compared to the original protein data. Here, we address the question whether it can still be used to reconstruct the regulatory structure and dynamics of the gap gene system in a manner which is consistent with previous efforts based on modelling, as well as genetic and molecular approaches to study gap gene regulation.
Before we present our modelling results, we provide a quantitative characterisation of our mRNA data, and compare them to the gap gene protein expression data described in [43]. Table 1 shows the number of embryos on which our mRNA data are based. Figure 3 illustrates the quality and resolution of embryo images underlying the two data sets. It shows time series of mRNA expression patterns for the trunk gap genes hb, gt, Kr, and kni produced with a colorimetric (enzymatic) in situ hybridisation protocol (columns 1–4 on the left), in comparison to protein expression data for Gt, Kni and the pair-rule protein Even-skipped (Eve) from the FlyEx database (http://urchin.spbcas.ru/flyex; [69], [70]; column 5 to the right). We used images such as the ones shown in columns 1–4 of Figure 3 to quantify the position of gap gene domain boundaries across space and time. This was done by applying the image-processing pipeline as described in Materials and Methods. In brief, we used texture-based image segmentation to create whole-embryo masks, which were used to rotate and crop embryo images. The developmental stage of each embryo was determined based on numbers of nuclei and membrane morphology as previously described [44] (see also Materials and Methods). We then extracted raw profiles of gene expression within a 10% strip along the embryo's lateral midline, and we determined the position of expression domain boundaries. The results of this analysis are shown in Table 2 and Figure 4. A detailed description of gap gene mRNA patterns can be found in Text S2.
As mentioned in the Introduction, we are mainly interested in the dynamics of gap domain boundary positions across space and time. Figure 4 compares those dynamics between mRNA and protein data. It is evident that mRNA expression patterns resemble those of proteins closely. The relative position of all domains with respect to each other is preserved at all time points. Anterior shifts in domain positions over time also mirror each other between mRNA and protein data: while mRNA domains are always more anterior than protein domains at equivalent stages (as reported in [46], Figure 3), the extent of the domain shifts is similar between the two data sets (see Table S1). Finally, domain sizes are similar as well (see Table S2), although protein domains of the gap genes are slightly larger than those of their transcripts.
However, there are also notable differences between mRNA and protein data. First, all mRNA patterns arise earlier than those of their corresponding proteins: mRNA expression of all gap genes is initiated before C13 ([40], and references therein), and gap mRNA domains are well established at early cycle 13. In contrast, protein levels have only just begun to be detectable, and increase rapidly, during that stage [43]. This is due to the delay caused by mRNA processing, nuclear export, and translation. Similarly, there is an evident lag between shifting positions of gap domains (see also Figure 3 in [46]), again indicating a significant delay between the dynamics of mRNA and protein patterns. In summary, while overall expression dynamics are similar between mRNA and protein, both the timing of expression and the absolute positions of gap domain boundaries differ between the two data sets.
The main aim of this study is to show that the gene circuit method—originally developed for protein expression data—can be adapted to work successfully with expression profiles derived from mRNA. To achieve this, gene circuit models were fit to mRNA expression patterns derived from our data set of boundary positions for hb, gt, Kr, and kni. External inputs to the model—regulatory contributions from maternal gradients Bcd and Cad, as well as terminal gap genes tll and hkb, which are not themselves regulated by trunk gap genes—were calculated based on protein data as described in Materials and Methods.
Our models run from early cleavage cycle 13 (C13), when gap proteins start to accumulate, to the end of C14A, when gastrulation starts. They span the trunk region of the embryo from 35 to 87% A–P position. This region is located slightly more anteriorly and is somewhat smaller than that used previously for protein models (35–92% A–P position [46], [48]–[50]). However, it covers the equivalent set of gap gene expression patterns, since mRNA domains are located more anteriorly and are slightly less wide than their corresponding protein domains (see Figure 4, and Tables 2 and S2). Fitting gap genes to the original range of 35–92% A–P position yielded equivalent results to those reported below (data not shown).
The quality of a gene circuit depends crucially on the relative position and dynamics of gap domain boundaries, while levels of expression are of secondary importance (see Introduction). However, there are two problems with the loss of information on gap gene product concentrations due to the way our data were processed. The first is that, after data normalisation, boundaries in our data set appear suddenly, without gradual build-up of gene product levels (see, for example, Figure 2O). This is clearly unrealistic, and leads to problems with the numerical stability of our gene circuit models, since very short time scales (and hence, very high production rates) will be favoured by the model-fitting procedure (data not shown). We have addressed this problem by re-scaling the data over time using a second-degree spline with a peak of expression at early C14A, capturing the gradual accumulation of mRNA during C13 and early C14A, as well as its degradation during late C14A ([40], [43], and our unpublished data). The second problem concerns the relative timing and intensity of expression in different gap gene domains. While central gap domains arise early, and show intense staining levels by C14A, more terminal ones arise later, and show less intense staining [43]. For instance, the anterior mRNA domain of hb arises at C9 or C10 from the maternal Hb gradient, while the posterior hb domain only appears at C13, with much lower intensity of expression [71]–[74]. We have addressed this problem by re-scaling the data along the A–P axis, using a scaling factor of 1.0 in the middle of the embryo (at 50% A–P position) that linearly decreases to 0.5 at the poles (0 and 100% A–P position). Note that this scaling step is justified by our prior knowledge of the gap gene system, but is not strictly required to generate the results we present below (see Discussion).
Finally, data were multiplied by a factor of 200, and collected into 50 bins (at C13) and 100 bins (at C14A) along the A–P axis to facilitate comparison with models based on protein data (see Materials and Methods and Figure S1 for details).
We used both ordinary (OLS) and weighted least-squares (WLS) fits to our full mRNA data set (see Materials and Methods). On one hand, the WLS method has the advantage of penalising ectopic expression outside the observed gap domains, and has been shown to be more effective than OLS for protein-based circuits [49]. On the other, it requires us to calculate the variance of each data point in our data set, which is not an obvious task given our data-processing methods. For this reason, we had to estimate weight values for the WLS cost function based on the observation that variance is higher in regions with high levels of gene expression than in those that show low, or now expression at all (see Materials and Methods for details).
We performed 150 fitting runs each with OLS and WLS cost functions respectively. Because of potential artefacts caused by overfitting, one cannot simply select those circuits with the lowest residual scores for further analysis. Instead, we inspected all solutions visually to detect obvious patterning defects (observed defects are described in Text S1). Only solutions without defects were selected. Using the WLS cost function resulted in a much higher fraction of runs with good fits to the mRNA expression data compared to OLS. For OLS solutions, only 10/150 (6.7%) were suitable for further analysis (Figure S2). Their residual errors—as measured by the RMS score (defined in Materials and Methods)—range from 16.6 to 20.9, with a median value of 17.453. For WLS solutions, 52/150 (34.7%) could be retained (Figure 5A). As is usual for WLS fitting [49], they show slightly higher RMS scores—between 20.0 and 21.7, with a median value of 21.052.
Overall, the selected gene circuits are extremely similar for both OLS and WLS (compare Figure 5A to Figure S3). Both show correct relative timing and positioning of gap domains, and reproduce the positional shifts of posterior domains towards the anterior. However, there are slight inaccuracies concerning the appearance and placement of domain boundaries. The formation of all gap domains, but particularly the posterior hb domain, is slightly delayed (Figure 5A, T1), and several boundaries are offset by 1–3 nuclei at specific points in time compared to the data (see, for example, the anterior boundary of kni at T1, and its posterior boundary at T3 in Figure 5A). These slight defects result in patterns that reproduce the data less faithfully than those obtained by WLS fits to protein data (Figure 5B; models from [49]). This is reflected in the corresponding RMS scores: 10.5 for protein data versus 16.6, the best score from our mRNA fits. In addition, our mRNA-based circuits show increased variability in model output between solutions compared to protein-based models (Figure 5, Figure S3).
Apart from these minor differences, however, there is significant agreement between all three sets of gene circuits. This similarity in expression patterns is reflected in the parameter values of our models. Distributions of estimated parameter values for regulatory weights of WLS-mRNA and WLS-protein solutions are shown as scatter plots in Figure 6A. Corresponding genetic interconnectivity matrices are compared in Figure 6B. It is clear from inspecting these matrices that a large majority of interactions are qualitatively the same in mRNA- and protein-based circuits (although increased variability in model output is reflected in increased variability of estimated parameter values for mRNA circuits). Repression corresponds to repression, and activation to activation, for a majority of circuits in both cases. Only six gap-gap cross-interactions are predicted to be different (emphasised by black frames in Figure 6B, and discussed in detail in Box 1), while all external inputs (bcd, cad, tll and hkb) interact with gap genes in the same manner for both data sets.
We analysed the differences between mRNA- and protein-based gene circuits in detail. Surprisingly, one of the differences is an improvement: our mRNA circuits predict repression of hb by Kr, which is in accordance with the experimental literature [77], [95]–[97]. This interaction had been missing in protein-based circuits. On the other hand, the two differences affecting the auto-regulatory terms of Kr and gt are of no functional significance, since auto-regulation is dispensable for correct gap gene expression [50].
The final three differences require more explanation: mRNA circuits predict Kr to inhibit kni, Kni to not interact with gt, and Gt to repress hb, while protein-based circuits predict no interaction for Kr on kni, and activation in the other two cases (Figure 6B). Of the mRNA interactions, only the model outcome of no interaction between Kni and gt is supported by experimental evidence, as the anterior boundary of the posterior gt domain is not affected in kni mutants [98]–[100]. With respect to the other two interactions, we know the effect of Kr on kni is an indirect one via the de-repression of gt. Also, the posterior hb domain appears normal in embryos over-expressing or mutant for gt [99], [101], [102]. In accordance with this, both mRNA- and protein-based circuits predict that the effect of Kr on kni, Kni on gt, as well as that of Gt on hb, is extremely weak (Figure 6A, and [48], [49]). Note that about 40% of mRNA-based circuits show an inhibitory effect of Kni on gt (Figure 6B). This suggests that it does not matter too much whether these interactions are activating or repressing as long as their strength is severely limited.
Interactions between overlapping gap domains are involved in regulating the dynamic anterior shifts of posterior expression boundaries over time. Jaeger et al. [46] suggested that the observed shifts are due to asymmetric repressive interactions, where the neighbouring domain at the posterior represses its adjacent domain at the anterior, but not vice versa. Therefore, repression of posterior by anterior neighbours (gt by Kni, and hb by Gt) poses a problem for this postulated mechanism. Still, our mRNA circuits exhibit correct shifting dynamics of posterior gap gene expression (Figure 5A). Graphical analysis of our mRNA-based models (Figure 6C) reveals the following relaxed condition for the shift mechanism: instead of complete absence of anterior-to-posterior repression it only requires an asymmetry in repressive strength. To be more precise, we define the net effect of an interaction as its regulatory weight multiplied by the concentration of the regulator integrated over the relevant region of the embryo (Figure 6C). Anterior shifts occur as long as the net effect of the anterior-to-posterior interaction is less repressive than the reciprocal posterior-to-anterior effect. This implies that the anterior-to-posterior interaction can be either activating or repressing, as long as it does not cancel out repression by the posterior neighbour. This allows the posterior repressor to become expressed in the posterior region of its anterior neighbour's domain, down-regulating the latter, and thus leading to the observed anterior shift of the latter's boundary.
In summary, none of the differences between mRNA- and protein-based circuits shown in Figure 6B are inconsistent with regulatory mechanisms for gap gene regulation postulated previously [46], [48]–[50] (see Box 1). The network recovered from both data sets is essentially equivalent. Our models—just as those obtained with protein data—predict that gap genes are regulated by (1) broad activation by maternal gradients, (2) auto-activation terms, which are non-essential, but serve to maintain sharply defined domain boundaries [50], (3) strong mutual repression between non-overlapping gap genes (hb and kni, Kr and gt; we call this mechanism ‘alternating cushions’ for reasons explained in [35]), (4) weak repression between overlapping gap domains showing posterior dominance which leads to dynamic anterior shifts in domain position, and (5) additional repression of gap genes at the posterior pole by terminal gap genes tll and hkb (Figure 6D; see [35] for review).
Sets of parameter estimates based on reduced-quality mRNA data show increased variability compared to protein-based circuits—both in terms of the distribution of parameter values (Figure 6A), and the regulatory categories they fall into (Figure 6B). Nevertheless, we have shown that we can recover consistent regulatory mechanisms from these estimates, if we consider the consensus network structure, that is, those regulatory categories into which a majority of the estimated parameter values fall into.
In this section, we examine if our parameter estimates are also determinable in the statistical sense defined by Ashyraliyev et al. [49], [64]. This is also known as practical parameter determinability analysis. It is achieved by assuming that our fitting problem has a single optimal solution, which we call the ‘true’ solution. We can then calculate confidence intervals around each one of our parameter estimates. These intervals determine a range of parameter values, which include—with a given probability of 95%—the true solution to the problem. A parameter is determinable, if its confidence interval lies entirely within one of our three regulatory categories: repression (parameter value less than −0.005), no interaction (between −0.005 and 0.005), or activation (greater than 0.005). It is weakly determinable if its interval intersects two of these categories, but excludes the third. There are two different ways to calculate these confidence intervals: dependent intervals tend to underestimate the extent of the confidence region, while independent intervals have the tendency to overestimate it (see Materials and Methods for details). As in previous studies [49], [64], we use independent intervals to estimate the determinability of a parameter.
Under specific conditions, gap gene circuits fit to protein data can yield parameter estimates, which are very well determined (Table 3, top row). Specifically, estimates from WLS fits are much more determinable than those obtained by OLS, if parameter values for diffusion rates (D), threshold parameters (h), and certain regulatory weights (the effect of Hkb on gt, Kr, and kni) are fixed during optimisation [49]. In contrast, determinability of parameter estimates is poor in both OLS and WLS fits to mRNA data (Table 3, rows 5 and 6). This loss of determinability could be due to three reasons: (1) the use of mRNA instead of protein data; (2) our data processing method (and approximation of variances for WLS); and (3) the smaller embryo region covered by mRNA-based circuits.
To distinguish between these three possibilities, we performed several series of optimisation runs. First, we used the original protein data set with approximated variances that only depend on expression level (as described for our mRNA data in Materials and Methods). This yielded a level of parameter determinability, which was somewhat worse, but still comparable to the original protein data set (Table 3, row 2). We then approximated expression boundaries in the integrated protein data set by our spline-based method, to emulate domain and boundary shapes equivalent to the mRNA-based data. Again, this yielded decreased but still reasonable parameter determinability (Table 3, row 3). Based on this, we conclude that neither the use of approximated weights, nor the use of approximated boundaries can account for all the loss of determinability observed in mRNA-based circuits. Finally, we performed fits to mRNA data on an extended region of 58 nuclei, equivalent to the region used for protein-based gene circuits. As for other mRNA-based solutions, determinability was poor, significantly decreased compared to any of the protein fits (Table 3, row 4). Taken together, the above evidence indicates that the most relevant factor for the loss of statistical determinability is the use of mRNA- instead of protein-based expression data.
In the previous sections, we have established that it is possible to reverse-engineer a developmental gene regulatory network from mRNA expression profiles. Even though these mRNA-based networks lack statistical determinability, the observation that we consistently recover a qualitatively equivalent network generates a certain confidence in the methodology. Next, we wanted to define the minimal requirements for gap gene circuits in terms of experimental data. As mentioned in the Introduction, the most time-consuming and hence limiting step involved in this method is the acquisition and processing of the quantitative experimental data required for model fitting. Thus, we would like to reduce and simplify our data as much as possible. To this end, we have performed a series of model fitting runs with reduced data sets (either by decreasing the number of boundaries, or the number of time classes), or with artificially generated data for external inputs (to simulate the case where quantitative protein data on maternal gradients and terminal gap genes would be unavailable). Table 4 presents an overview of the distinct data sets used for model fitting, while the analysis of the resulting gene circuits is summarised in Figure 7 (based on interconnectivity matrices shown in Text S3).
Our results show that a smaller fraction of solutions obtained by fitting to reduced data sets are usable for analysis. This applies to both reduction of the number of boundaries and time points present in a data set (Figure 7, 2nd column). In contrast, the optimum and median RMS scores of these fits do not show any clear trend, although RMS scores slightly increase as boundaries are removed (Figure 7, 4th column). Using artificial external inputs did not significantly affect the number of usable solutions or the RMS score. We conclude that our approach is surprisingly robust to the reduction in quality of our data set. Only in one case—reduction of the number of boundaries to 20%—did we completely fail to get any solutions suitable for further analysis.
While a good fit to the data is important, it is even more crucial that gene circuits represent consistent gene network structures and regulatory mechanisms. For this reason, we examined whether models fit to reduced or artificial data incorporate the five basic mechanisms of gap gene expression described in the previous section.
First, we looked at mutual repression of non-overlapping gap genes (alternating cushions). This mechanism is recovered extremely robustly: it is present in all selected solutions (Figure 7, blue column). In agreement with previous analyses [50], we conclude that it lies at the core of gap gene regulation.
With respect to auto-activation, we observe that all mRNA-based circuits behaved similarly, with the exception of the gene circuits from fits with artificial maternal gradients (Table 4, and Figure 7, green columns). In general, hb, gt and kni were auto-activating, while Kr was not (or not consistently, at least). In addition fits with artificial maternal gradients show reduced presence of gt auto-activation. The variability in these results can be explained by the fact that auto-activation is not essential for positioning gap domains in gene circuit models [50].
To analyse the presence or absence of the domain shift mechanism [46], we examined net effects of repression between overlapping gap domains as described in the previous section (Figure 7, red columns; see also Figure 6C). We observed that two domain interactions were consistently present in our models: the net repression of Kr by Kni, and that of kni by Gt. However, net repression of gt by Hb was only recovered in about 10–50% of solutions. In those circuits that do not show this effect, Gt represses hb in the anterior to reproduce the peak of hb expression in the middle of the embryo. We have shown elsewhere that this is a modelling artefact [48]. Repression of hb in the posterior is overcome in these circuits by strong activation of hb by Tll (data not shown), an interaction which is likely to be indirect, and therefore not supported by experimental evidence [75].