Efficient reverse-engineering of a developmental gene regulatory network
Unlike the interactions described above, mutual repression between the overlapping domains of hb and Kr does not contribute to gap domain shifts [46]–[48]. If we consider the net effect of these reciprocal interactions, hb is repressed by Kr in the majority of cases (Figure 7, orange column). In contrast, most protein-based circuits show net repression of Kr by Hb. We explain this as follows: the border between hb and Kr is maintained at the same A–P position over time (Figure 4; see also [43], [46]). There are several regulatory mechanisms that can accomplish this. Protein-based circuits exhibit auto-activation on both sides plus a minimal net repressive effect of Hb on Kr, while mRNA-based circuits favour auto-activation of hb only, plus stabilisation of the hb/Kr border by a slight net repression of Kr on hb (Figure 7, green and orange columns). The experimental evidence on which of these alternative mechanisms applies in the embryo remains inconclusive: while the potential presence of Kr auto-activation [76] favours the first, the absence of any defect in the posterior hb boundary in Kr mutants supports the second mechanism [77], [78].
Next, we looked at activation of gap genes by the maternal factors Bcd and Cad. This mechanism is recovered very robustly (Figure 7, purple column), in agreement with previous studies [46], [48]–[50]. Note that the summary bar graphs of Figure 7 (purple column) omit the effect of Cad on hb and Kr, as these interactions are an artifactual aspect of the model ([48], and references therein).
Finally, the terminal gap genes Tll and Hkb repress gap genes at the posterior pole of the embryo. In general, we recover this mechanism robustly (Figure 7, brown column). However, if Tll and Hkb protein gradients are substituted by mRNA profiles (see Table 4), about 20% of the solutions fail to show repression of Tll on gt. In these circuits, repression of gt by Hb is sufficient for the retraction of its posterior domain from the pole (data not shown). Interestingly, recovery of this interaction is rescued in those fits in which both maternal and terminal external inputs were replaced by ‘artificial’ patterns (see Figure 7).
In summary, we have managed to recover surprisingly correct and accurate gap gene regulatory mechanisms even with data sets of severely reduced coverage and/or quality. This does not imply that we can reconstitute specific gene networks with arbitrary data. On the contrary, successful network inference requires very specific conditions for the data used in model fitting. We will revisit this important point in the Discussion.
Earlier studies using reverse-engineering approaches were based on quantitative protein expression data [46], [48]–[50]. These protein data capture the precise shape of each expression domain, as well as differences and changes in relative protein levels both between domains of the same gene, and specific gap domains over time. In this paper, we have demonstrated that we can correctly infer the regulatory structure and dynamics of the gap gene network in Drosophila melanogaster using spatial mRNA-based expression data of reduced quality and coverage. Our mRNA data only capture the dynamics of boundary positions, and even that at a reduced spatial resolution. Furthermore, our simplified data processing pipeline leads to a loss of information on relative expression levels, which need to be approximated based on straightforward biological assumptions on transcriptional regulation and gene expression (gradual increase and decrease of levels over time).
On the other hand, this type of mRNA data can be acquired and processed within a fraction of the time and effort it takes to obtain high-quality protein-based expression patterns: instead of years, it now only takes a few months to establish quantitative expression profiles for a complete set of genes involved in a developmental process. Our approach avoids having to raise antibodies against regulators, which is both technically challenging and expensive. It uses a robust colorimetric (enzymatic) staining protocol instead of fluorescence. It avoids laborious scanning of embryo images using confocal microscopy (which is not available everywhere). In addition, we show that much fewer embryos need to be processed and stained, and fewer developmental stages need to be represented in the data than previously thought. And finally, our approach simplifies data processing, further reducing the effort required for data quantification.
Does the fact that we still recover the correct gap gene network using mRNA data imply that our method lacks specificity? Would it infer the same network with any kind of data? The following evidence demonstrates that the answer to these questions is a clear no, and suggests minimal conditions for the expression data that must be met for inference to be specific and consistent.
First of all, gap gene circuit models fail to correctly predict gap gene expression in the head region of the embryo (anterior of 35% A–P position) [45], [48]. In this region, additional regulators (the head gap genes which are not included in our models; see [35] for review) are required for correct regulation and expression. This implies that the model-fitting procedure is specific: it fails when relevant regulatory factors or mechanisms are missing.
In contrast, we have mentioned earlier that accurate measurements of absolute expression levels do not seem to be crucial for correct network inference (see Introduction). Our results suggest that we can recover the structure of the gap gene network even if relative levels are only approximated in the data. This confirms that the most important expression feature for network inference is the dynamic positioning of expression boundaries.
This is further corroborated by the following: early attempts at reverse-engineering the gap gene system using gene circuits exhibited various patterning defects and showed inconsistent network structures between independent model fits [45]. In contrast, we are able to infer correct network structure using data sets with spatial and temporal resolution similar to that used in [45] (Figure 7). Moreover, our data sets approximate gap protein domain positions with mRNA data (Figure 4). This suggests that the failure of inference in the earlier study was not caused by a lack of accuracy or resolution, but rather by qualitative errors regarding the relative placement of domains in the data used for model fitting. As mentioned in the Introduction, the abdominal kni domain was placed too far anterior in those data, leading to an exaggerated overlap with the central Kr domain, and an artefactual gap between the domains of kni and hb in the posterior of the embryo (see Figure 1 in [45]). Similar problems affected an early study of the eve gene: in the data used for fitting, the position of the fifth expression stripe was shifted compared to the domains of its gap gene regulators, which led to inaccurate prediction of the regulatory mechanism underlying its expression [39].
What differs between these early attempts and more recent reverse-engineering studies is that both protein- [43] and mRNA-based data used in the latter capture the relative arrangement and timing of gap domains correctly. Both data sets show qualitatively equivalent expression patterns (see Results, Figures 3 and 4). The order in which gap domains are arranged along the A–P axis, as well as the order in which they appear relative to one another are the same for mRNA and protein data. Evidently, this temporal and spatial order is determined by the specific regulatory structure of the network (given suitable kinetic parameters). This, in turn, allows us to robustly recover the specific regulatory structure of the gap gene network as long as our data capture their particular kind of spatio-temporal dynamics.
Obviously, there are limits to the amount of inaccuracy the data can contain with regard to absolute positions and variability of the patterns. This is shown by the fact that we were unable to recover any usable gap gene circuits with data sets reduced to 20% of the original number of measured domain boundaries (Figure 7). These data sets only contain about 1–3 embryos per boundary per time class, while the smallest successful data set (40% data) has about 1–7 embryos on average.
In general, the number of embryos that need to be processed will depend on the natural variability of the patterns and the quality of the experimental protocols: noisy expression data require larger sample sizes. Similarly, temporal resolution will depend on the time scale of patterning dynamics. In our case, only 3 time classes (each about 25 min apart) were sufficient to recover a correct regulatory network. This is probably due to the fact that gap domains shift and develop smoothly over time (see Figure 4). Systems with more abrupt or uneven changes in expression patterns will require a higher temporal resolution.
Our investigation of minimal data requirements for model fitting is of an empirical nature. It should be corroborated and extended in the future by more systematic and rigorous approaches based on methods for optimal experimental design (OED, reviewed in [79]). OED uses algorithms for global optimisation and calculation of confidence intervals similar to those used here to predict which measurements in a data set are most relevant to accurately estimate parameter values (see, for example, [80], [81]). This could be used for a more rational design of data sets used for reverse-engineering by guiding the choice of observables and time points that are most informative to infer network structure and dynamics. However, applying OED to real-world, complex, non-linear systems remains challenging, and has only be achieved in exceptional cases (see [82]). Therefore, it is beyond the aim and scope of our current study.
In summary, our results demonstrate that the amount of data required for reverse engineering is much lower than previously thought. The necessary data sets can be acquired and processed by a single researcher within a few months, without the need for expensive equipment. The precise number of embryos to be processed, and the required temporal resolution need to be adapted according to the features of the system under study, the main condition being that the data capture the relative timing and spatial arrangement of expression domains correctly.
The results reported in this study—together with earlier evidence from protein-based circuits [49]—indicate that weighted least squares (WLS) fits perform much better than ordinary least squares (OLS). In principle, however, WLS fits require the accurate measurement of variances in expression levels to calculate the weights for the sum of squares. This is not possible in our current data quantification framework. Moreover, in general it increases the amount of work required for data acquisition and processing considerably. First, methods that measure the relative level of expression accurately are technically more challenging and require more work than those presented here (see previous section). And second, enough individual expression patterns need to be quantified for the measured variance to be reliable.
We avoid this complication by using approximated weights for WLS, which are simply proportional to the expression level of a gene. This assumption is mainly based on methodological reasoning (although it is also supported by experimental evidence on protein expression patterns [40], [43]). It is crucial to avoid small ectopic expression domains in non-expressing regions that can exert significant regulatory effects if interaction weights are sufficiently large. Several such ectopic domains have been observed in fits using OLS (our data and [47], [49], [51]). WLS fits effectively suppress these modelling artefacts, as long as the penalty for ectopic expression imposed by small variances in non-expressing regions are sufficiently high. This suggests that our choice of approximated variances is justified for practical reasons, since it emphasises the importance of boundary positions during the fitting process.
Although we do recover the same network from our mRNA-based models as that predicted by protein-based circuits, there is much increased variability in the estimated parameter values (Figure 6A,B). This is also reflected in the loss of parameter determinability we observe in our results (Table 3). In practice, both of these phenomena imply that our mRNA-based models do not predict one, but rather a small set of possible network structures, while protein-based models predicted a specific, single network (see Figure 6B). This is not surprising, since the quality of a gene circuit model reflects the quality of the data it was fit to. But of course it is a problem, although, we would argue, not a fundamental one. It can easily be addressed by combining the reverse-engineering approach with experimental (genetic and molecular) verification.
Let us illustrate this with an example: there are several alternative network variants that occur in our circuits. They are all minor in that they differ from each other in one or two interactions at most. One of these variants occurs in models in which the posterior domain of hb arises due to strong and direct activation by Tll, rather than the absence of repression by Gt, its immediate anterior neighbour (see Results). This alternative mechanism is a plausible explanation for the observed expression dynamics. However, it is not compatible with experimental evidence, which enabled us to classify it as an artefact of the model.
Note that both the experimental evidence (reviewed in [35]) and predictions based on gene circuits [48], [49] contain such ambiguities. Fortunately, these unresolved (and potentially not resolvable) cases have only a limited overlap, since genetic and reverse-engineering approaches are complementary to each other: one inferring regulatory interactions from mutants, the other from wild-type expression patterns [48]. In other words, interactions which are not clearly supported by experimental evidence are mostly unambiguous in our gene circuits, while others that are ambiguous in our models (as, for example, the repression of hb by Gt), are clearly supported (or excluded) based on experimental evidence.
But what about systems which have been less well studied than Drosophila? In such systems, we do not have a comprehensive experimental literature to compare our results to. Instead, model predictions will have to be tested using genetic approaches such as mutant analysis, over-expression assays, or gene knock-down by RNA interference (RNAi). Again, the reverse-engineering method is most powerful when used in conjunction with complementary experimental approaches.
By minimising the amount of quantitative data required for reverse-engineering a developmental gene regulatory network, we have removed a major bottleneck for applying the method more widely. Still, this method is unlikely to be scalable to systems that are orders of magnitude larger than the one studied here. Microscopy and image acquisition remain labour-intensive, and our quantification pipeline still requires a series of manual interventions, such as positioning the splines that are fit to expression boundaries, or time classification of embryos. It remains a major challenge to fully automate these steps. Therefore, the effort required to quantify hundreds or thousands of spatial gene expression patterns is still considerable, even if robust and fast methods are used (see, for example, [83], [84]). Moreover, global non-linear optimisation is computationally intensive, and may not yield unique solutions in large regulatory systems.
On the other hand, many pattern-forming networks are similar in complexity and nature to the gap gene system. One example in Drosophila is the dorso-ventral patterning system in the early embryo [85]–[87]. Other developmental networks occur in the context of cellularised tissues and often involve more than just one spatial dimension. It is straightforward to extend the gene circuit method to such systems. Complicating factors, such as post-transcriptional regulation, cell-to-cell signalling, or tissue movements and growth, can readily be accommodated in the gene circuit modelling formalism [38]. There are several examples of tractable cellular patterning systems in Drosophila: mesoderm and heart development [88], [89], morphogen-based patterning in wing imaginal discs [87], [90], or the thoracic bristle patterning system [91], [92]. Examples beyond Drosophila include vulval induction in the roundworm Caenorhabditis elegans (see, for example, [93]), or the dorso-ventral patterning system of the vertebrate neural tube [94].
Generalisation of the method is further facilitated by the fact that the approximations we have used are based on straightforward assumptions that do not require any in-depth understanding of the system under study. Our approximation of relative gene expression levels assumes smooth increase and decrease of gene product levels. The scaling of expression levels from the centre to the terminal ends of the embryo is a reasonable assumption in the case of Drosophila, but may be omitted in any other system as it is not strictly required for inferring the correct network structure. Approximated variances simply assume low levels of variability in non-expressing regions. Artificial external inputs were based on qualitative inspection of maternal gradients and expression of terminal gap genes. All of these assumptions can be applied to systems other than the gap gene network. The only potentially problematic assumption concerns the use of mRNA expression profiles. Obviously, one should have some qualitative evidence that mRNA patterns indeed resemble protein profiles before relying on such data.
Application of reverse engineering to a large number of developmental systems in different organisms would allow us to investigate the dynamics of pattern-forming gene networks in a quantitative and systematic manner [26]. Our simplification of the method indicates that this is achievable within a reasonable amount of effort. The potential benefits of such a research programme are significant. Only through quantitative investigation of specific, experimentally accessible, gene networks will we be able to better understand the principles that govern development and pattern formation.
A particularly interesting application of the gene circuit method is the comparative analysis of homologous gene networks across different species [26]. Such a comparative analysis allows us to identify conserved and divergent regulatory interactions in an evolving network. Moreover, changes in regulatory mechanisms can be mapped—rigorously and unambiguously—to observable differences in gene expression between species. We are currently performing such an analysis between gap genes in Drosophila and other dipteran species. It will allow us, for the first time, to study the evolution of a real developmental gene regulatory network in a detailed and quantitative manner. This is a crucial step towards a more general investigation into the causal relation between evolution at the genotypic and the phenotypic level.