Efficient reverse-engineering of a developmental gene regulatory network

Understanding the complex regulatory networks underlying development and evolution of multi-cellular organisms is a major problem in biology. Computational models can be used as tools to extract the regulatory structure and dynamics of such networks from gene expression data. This approach is called reverse engineering. It has been successfully applied to many gene networks in various biological systems. However, to reconstitute the structure and non-linear dynamics of a developmental gene network in its spatial context remains a considerable challenge. Here, we address this challenge using a case study: the gap gene network involved in segment determination during early development of Drosophila melanogaster. A major problem for reverse-engineering pattern-forming networks is the significant amount of time and effort required to acquire and quantify spatial gene expression data. We have developed a simplified data processing pipeline that considerably increases the throughput of the method, but results in data of reduced accuracy compared to those previously used for gap gene network inference. We demonstrate that we can infer the correct network structure using our reduced data set, and investigate minimal data requirements for successful reverse engineering. Our results show that timing and position of expression domain boundaries are the crucial features for determining regulatory network structure from data, while it is less important to precisely measure expression levels. Based on this, we define minimal data requirements for gap gene network inference. Our results demonstrate the feasibility of reverse-engineering with much reduced experimental effort. This enables more widespread use of the method in different developmental contexts and organisms. Such systematic application of data-driven models to real-world networks has enormous potential.

Only the quantitative investigation of a large number of developmental gene regulatory networks will allow us to discover whether there are rules or regularities governing development and evolution of complex multi-cellular organisms.

Elucidating the regulatory structure and dynamics of gene networks is a major objective in biology. The inference of regulatory networks from gene expression data is known as reverse engineering [1]–[7]. It is being widely and successfully applied, from microbes to animals (see, for example, [8]–[16]). Many reverse engineering studies aim to determine regulatory structure from large-scale perturbation- or time-series data based on microarray or transcriptome-sequencing technology (reviewed in [5], [17]). This approach has two significant limitations: first, spatial information on gene expression is lost, since homogenised tissue samples or disaggregated cells are studied. And second, most resulting models are of a static and probabilistic nature, which cannot be used to investigate network dynamics (for example [18]–[20]). If dynamical models are used, they are often linear (for example [21]–[23]). Network inference using complex non-linear dynamical models is deemed a considerable technical challenge [6], [17], [24]. However, there are many important biological questions that absolutely require consideration of non-linear and spatial aspects of a system. Here we discuss such a case, and show that reverse engineering can be used for its study with a reasonable amount of experimental and computational effort.

Our research focuses on how developmental gene regulatory networks produce spatial patterns in multi-cellular organisms, and how these patterns evolve through changes in the underlying structure of the network [25], [26]. In this context, reverse engineering is implemented by fitting non-linear systems of differential equations to quantitative, spatial gene expression data (reviewed in [7]). There are many network modelling formalisms [27]–[29], and a number of powerful global non-linear optimisation methods [7], [30], [31], which are suitable for this task.

So far, only a small number of developmental systems have been reverse-engineered using dynamical models (see, for example, [32]–[34]). One of those is the (trunk) gap gene network of the vinegar fly Drosophila melanogaster (reviewed in [35]). This regulatory network consists of four genes—hunchback (hb), Krüppel (Kr), giant (gt) and knirps (kni)—which all encode transcription factors. Gap genes are involved in establishing the segmented body plan of the animal. They are active during a very early period of Drosophila development, called the blastoderm stage, which occurs before the onset of gastrulation. At this stage, the embryo consists of a multi-nucleate syncytium allowing transcription factors to diffuse through the tissue. Gap genes are expressed in broad, overlapping domains along the embryo's major, or antero-posterior (A–P) axis. They are regulated by long-range gradients of transcription factors encoded by maternal co-ordinate genes bicoid (bcd), hunchback (hb), and caudal (cad), and are repressed by the terminal gap genes tailless (tll) and huckebein (hkb) in the pole regions of the embryo. Maternal co-ordinate and gap genes form the first two tiers of the segmentation gene hierarchy in Drosophila. Together they regulate pair-rule and segment-polarity genes, the latter forming a molecular pre-pattern that leads to morphological segmentation at later stages of development (see, [36], [37], for review).

The particular reverse-engineering approach we use to investigate the gap gene network is called the gene circuit method [1], [38], [39] (Figure 1). It uses mathematical models called gene circuits that represent the basic properties of the embryo and the transcriptional regulatory interactions underlying the network. Gene circuits are described in detail in Materials and Methods. Here we provide a brief overview of the model, which consists of a row of dividing nuclei (Figure 1, top left) each harbouring an identical version of the regulatory network. There are three processes that occur within and between nuclei: (1) regulated gene product synthesis, (2) Fickian gene product diffusion, and (3) linear gene product decay (Figure 1, top middle). Regulatory interactions that direct synthesis are represented by a genetic interconnectivity matrix: each regulatory weight in this matrix can represent activation, repression, or no interaction depending on whether it is positive, negative or (close to) zero (Figure 1, bottom panel, right). The interconnectivity matrix can also be displayed as a network diagram (Figure 1, bottom panel, left). Note that the values of regulatory weights are not set a priori. Instead, they are estimated by fitting the model to quantitative gene expression data. Those model solutions that fit the data well are analysed to characterize the regulatory structure and dynamics of the network (Figure 1, bottom). In this way, gene circuits act as tools to extract regulatory information from quantitative data.

Previous reverse-engineering studies of the gap gene network were based on quantitative expression data obtained by visualising the distribution of gap gene mRNA [40], or protein products [41]–[43] using fluorescent whole-mount in situ hybridisation, or antibody staining (immunofluorescence) respectively. Stained embryos were imaged using confocal laser-scanning microscopy, and the resulting expression profiles were quantified using a processing pipeline that includes image segmentation to identify nuclei, time classification, removal of non-specific background staining, data registration to remove embryo-to-embryo variability, and data integration (reviewed in [44]). It took years of effort by several researchers to establish the protein data set [41], [43]. mRNA data, on the other hand, were acquired by one of the authors of this paper in considerably less time [40]. However, these mRNA data remain incomplete in that they only cover a subset of gap genes (Kr, kni, and gt) during the earliest stages of expression.

These previous reverse-engineering studies have yielded many new insights into gap gene regulation, which would have been difficult to obtain by experimental approaches alone. An early pioneering study predicted a co-operative effect between maternal factors Bcd and Hb on the regulation of gap gene expression [45]. Later efforts uncovered a mechanism for the dynamic anterior shift of gap domains over time [46], [47], removed ambiguities in the published experimental evidence [48], [49], identified core mechanisms for gap gene regulation [48], [50], and explained the robustness of the system against variable levels of maternal inputs [47], [51]. Taken together, these studies clearly demonstrate the utility and feasibility of the approach: over the past two decades, reverse engineering has contributed significantly to our understanding of gap gene regulation.

It would be extremely interesting to apply the gene circuit method to other developmental systems. In our view, reverse engineering has tremendous potential for the study of gene regulatory networks in development and evolution. For instance, gene circuits could be used to reconstruct homologous developmental regulatory networks across a range of species, to compare their regulatory structure and dynamical behaviour [26]. This could be used to map which regulatory changes in a network correspond to which changes in gene expression during evolution. Alternatively, it would also be highly interesting to compare network structures and dynamics between different developmental processes.

Despite this potential, the application of dynamic, non-linear reverse-engineering approaches beyond gap genes in Drosophila has been very limited. The main reason for this, we suspect, is the following: collection of high-quality data sets—such as the spatio-temporal profiling of gap genes described above—is costly both in terms of time and resources. It is clearly the bottleneck of the approach. Protocols based on immunofluorescence require antibodies, which are difficult and expensive to obtain. Confocal microscopy is time-consuming and laborious, since a large number of embryo images need to be scanned. Moreover, while protocols for data acquisition and quantification work efficiently in Drosophila, their application to less well-established experimental models is not trivial. In particular, it is often difficult to adapt fluorescent staining protocols to non-model species.

Thus, in order to make the gene circuit method more widely applicable—and hence useful for the study of developmental gene regulatory networks—it is imperative that we simplify the method. We address an important question which applies to reverse-engineering approaches in general: how much, and what kind of data are required to successfully infer a gene regulatory network? Answering this question in the context of the gap genes will allow us to minimise the cost of data acquisition and processing. This, in turn, will decrease the barrier for applying reverse-engineering methodology to other developmental systems, many of which are similar in kind and complexity to the gap gene network.

The quality of a gene circuit model depends directly on the quality of the data it was fit to. What matters most in this regard is the timing and position of expression domain boundaries with respect to each other. The relative level of expression in each domain is less crucial. For instance, early gap gene circuit models did not capture the formation of the abdominal kni domain correctly (see Figure 2 in [45]). This was due to the incorrect relative position of this domain in the data resulting in a large gap between it and the posterior hb domain (see Figure 1, ibid.). This defect is no longer present in more recent models based on data with the abdominal kni domain positioned accurately while still only measuring relative levels of protein concentration [46], [48]–[50].

In this study, we present a simplified reverse-engineering protocol and apply it to a new, quantitative data set of gap gene mRNA expression in Drosophila. We demonstrate how mRNA expression data derived from a colorimetric (enzymatic) protocol for in situ hybridisation can be used to infer the regulatory structure and dynamics of the gap gene network. We compare our results with those obtained in previous studies based on protein expression data, and show that they predict equivalent regulatory mechanisms that are consistent with experimental evidence. In addition, we show that our simplified data set can be reduced even further while still yielding correct predictions. In this way, we define a set of minimal requirements for the successful inference of gap gene regulatory network structure and dynamics. These minimal requirements suggest that the adapted gene circuit method can be applied to a variety of developmental systems with a reasonable amount of effort. Such wider application of reverse-engineering methods will enable us to carry out systematic and comparative analyses of developmental gene regulatory networks.

cDNA clones were ordered from the Drosophila Genomics Resource Center (DGRC; dgrc.cgb.indiana.edu) and used to make riboprobes labelled with DIG and/or FITC. Wild-type blastoderm-stage Drosophila embryos were collected after 4 hrs of egg laying and stained with a colorimetric in situ hybridisation protocol adapted from [52] and [53]. In brief, fixed and dehydrated embryos were re-hydrated by washing 1x in PBT/methanol (embryos are allowed to sink before the solution is removed), 2x in PBT, and 1×5 min in PBT. Embryos were incubated in PBT containing 0.179U proteinase K for 1 hr on ice, then washed 2x in ice-cold PBT. Embryos were post-fixed for 25 min with 5% formaldehyde in PBT with mild shaking, then washed 1x followed by 2×5 min in PBT. Embryos were pre-hybridised by washing 1×10 min in equal volumes PBT and hybridisation buffer (50% formamide, 5xSSC, 5 µg/ml yeast tRNA, 100 mg/ml salmon-sperm DNA, 50 µg/ml heparin, 0.1% tween-20 in DEPC-treated water), 1×2 min in hybridisation buffer, and 1×1 hr in hybridisation buffer at 56°C. Hybridisation was carried out overnight: 0.5–1 ng/µl of probe(s) were added after heating at 80°C in a small amount of hybridisation buffer for 3 minutes. Post-hybridisation, the embryos were washed 1×15 min and 2×30 min in hybridisation buffer, 1×15 min in equal volumes of PBT and hybridisation buffer, and 4×15 min in PBT. Blocking steps were carried out using 5% heat-treated goat serum in PBT for 30 min, followed by incubation with anti-DIG or anti-FITC antibodies conjugated with alkaline phosphatase (Roche) at 1∶2000 in 5% heat-treated goat serum in PBT for 1 hr. Unbound antibody was removed with washes of 3x followed by washes of 4×15 min in PBT. To prepare for staining, embryos were washed 2×5 min in AP buffer (100 mM NaCl, 50 mM MgCl, 100 mM Tris pH 9.5, 0.1% tween-20).

Staining was carried out in the dark by the addition of AP buffer containing 0.1 mg/ml NBT and 0.05 mg/ml BCIP. Staining was stopped with 3x washes in PBT. For single staining (one probe), embryos were washed a further 3×10 min then counter-stained (see below). For double staining (two probes), alkaline phosphatase was inactivated by washing 1×1 min, then 1×10 min in glycine buffer (0.1 M glycine 0.1% tween pH 2) followed by 3×10 min in PBT. Blocking, antibody incubation and washing steps were carried out as described above. To prepare for staining, embryos were washed 2×5 min in Fast Red buffer (100 mM Tris pH 8.2, 0.1% tween-20). Staining was carried out in the dark by the addition of Fast Red solution (1 Fast Red tablet (Roche) dissolved in 2 ml Fast Red buffer). Staining was stopped with 3×1 min followed by 3×10 min washes in PBT. Nuclei were counter-stained by a 10-min incubation in PBT containing 0.3 µM DAPI, followed by washes of 3x followed by 3×10 min in PBT. Embryos were cleared through a series into 70% glycerol:PBT, of which 30 µl were mounted per slide. All washes were done on a nutator, except for those in proteinase K.

An overview of image acquisition and processing steps is shown in Figure 2. For each embryo, four images were acquired using a compound, wide-field, fluorescence microscope: (A) a bright-field image (Figure 2A), (B) a fluorescent image of the DAPI nuclear counter-stain (Figure 2B), (C) a differential interference contrast (DIC) image (Figure 2C), and (D) a DIC image of membrane morphology on the dorsal side of the embryo (Figure 2D). Images A–C were acquired using a 10x objective, image D using a 40x objective. Images A and B are focused on the surface, images C and D on the sagittal plane of the embryo. All images were taken at 8-bit accuracy, thus setting the range per RGB channel to [0,255].

Only laterally oriented embryos were selected for processing. Gene expression patterns were extracted from embryo images as follows. Binary masks covering the whole embryo are calculated using a sequence of image segmentation steps on the DIC image (Figure 2C). Intermediate steps are shown in Figure 2E/F: 1. the RGB image is converted to gray-scale, 2. gamma correction is applied to increase contrast, 3. the image is inverted, 4. Sobel edge detection is carried out, 5. dilation operations are applied to the resulting binary image, and 6. holes are filled [54], [55]. This results in a number of contiguous binary blobs in the mask image. All blobs touching the image border are removed, and only the largest blob is retained. Finally, a smooth whole-embryo mask is created by applying a Gaussian filter to the remaining blob. This mask and all raw images (A–C) were rotated and cropped as described in [56] such that the embryo's major, or antero-posterior (A–P) axis is horizontal. If necessary, embryo images were flipped manually to a canonical orientation such that anterior is to the left, and dorsal is up (Figure 2G–J).

To extract gene expression profiles from an embryo, a smooth cubic spline was generated with five equidistant knots through the main branch of the skeleton of the embryo mask (Figure 2K; [55], [57]). Along the spline, we extract average RGB-values over a band of 10% mask height (5% above and below the spline curve), resulting in raw profiles for each RGB channel of the bright-field image (Figure 2H) along the A–P axis (Figure 2K,L). Values for FastRed or NBT/BCIP staining were then calculated as follows: FastRed = green − red, and NBT/BCIP = red, where ‘red’ and ‘green’ refer to inverted RGB colour channels extracted from the bright-field image. Inspecting 1D-graphs of the resulting profiles (Figure 2M), boundaries for gene expression domains were extracted manually for each embryo: each boundary was labelled with a unique identification number (see Supplementary Material), and two points (x0, y0) and (x2, y2) were determined that indicate the beginning and end of the boundary, where staining levels approach background and maximum levels respectively. A middle, third control point (x1, y1) was automatically calculated from the other two points by taking the average for x, and locating the corresponding expression level y. Hence represent relative position along the A–P axis (in percent, where 0% is the anterior pole), and represent the relative intensity level of the staining. Points x0,2 and y0,2 were used as anchor points for cubic splines with fixed zero-derivatives at their end knots. Finally, splines were normalised such that the expression level at the starting point was 0, and the expression level at the end point was 1 (Figure 2N).

Integrated time-series of gene expression were prepared as follows. Embryos were staged into separate cleavage cycles (defined as the period between mitotic divisions n-1 and n; e.g. cycle 13 spans the time between mitoses 12 and 13) based on nuclear density and number of nuclei in images showing DAPI nuclear counter-stains (Figure 2B; time classification for early embryos was described previously in [40]). C14A was further subdivided into eight equally spaced time classes (T1–8) based on membrane morphology from high-resolution DIC images (Figure 2D; time classification for late embryos as described in [44]). Expression domain boundaries were grouped by gene, stage and boundary identification number (see above). Average boundary positions were determined by calculating separately the median start and end points for each group, which were then used for fitting a median-boundary spline as described for individual boundaries above. Finally, we combined different domain boundaries for each gene at each time class into an integrated, normalised expression profile along the A–P axis (Figure 2O). Spatial registration of domains was performed by checking the integrated expression profiles against double stained embryos (i.e. embryos stained for two gap genes, for instance hb and Kr) to verify the relative spatial order of gap gene expression domains (data not shown).

The following post-processing steps had to be applied to our data to make them suitable for model fitting and comparison (see Results, and Figure S1 of the online Supporting Information): (1) We collected our normalised, integrated mRNA expression data into 50 (C13) or 100 (C14A) bins to reflect the approximate number of nuclei along the A–P axis [46], [48], [49]. (2) We scaled the intensity of our expression data linearly along the A–P axis from 50% A–P position (mid-embryo; ×1.0) to both termini (poles of the embryo; ×0.5) to reflect the higher intensity of central versus more terminal gap domains [43]. For Drosophila expression data, this is a reasonable assumption, but not an essential requirement. Omitting this post-processing step resulted in qualitatively equivalent results (data not shown). (3) We also scaled our expression data along the time axis by a second-degree spline with a peak of expression during early C14A, to capture the gradual accumulation (during C12, C13 and early C14A) and degradation (during late C14A) of gap gene mRNA ([40], [43], and our unpublished data). Normalised boundaries were scaled to 0.1 at the onset of expression (early C13, t = 0.0 min), 1.0 at around T5 (t = 48.0 min), and 0.7 at gastrulation time (t = 71.1 min), which is the final time point. (4) We multiplied our mRNA expression data by a constant factor of 200. This makes the scale of both mRNA and protein data match as closely as possible, and therefore facilitates comparison to models obtained with Drosophila protein data [46]–[51].

The posterior Kr domain, which arises in late C14A, was removed from the data used for model fitting to avoid modelling artefacts. This domain is known to be under regulatory control of the terminal gap genes with additional inputs from the Forkhead (Fkh) transcription factor (not included in this study), and it does not participate in segment determination [58].

Fitting models using a weighted least squares (WLS) protocol (see below) requires a weight for each data point indicating its associated variation. As our mRNA expression data do not provide such information, weights were created from normalised, integrated mRNA expression data according to the formula: with being the normalised staining intensity and v the corresponding weight. This proportionality of variation with expression level reflects the fact that gap domains (showing high levels of expression) show more variation than those regions of the embryo in which a gene is not expressed [40], [43].

Gene circuits used for model fitting require external inputs based on data for maternal gradients (Bcd and Cad) as well as terminal gap genes (Tll and Hkb). Depending on the scenario we wanted to test (see Results), expression profiles for these inputs were taken from previously published quantitative protein data [43], [49], or were approximated as follows. Bcd: a time-independent (i.e. constant) anterior gradient was created by fitting an exponential curve to all available Bcd data across time and space. Cad: an artificial posterior morphogen gradient was created using 2D thin-plate splines [59] based on the following minimal set of features of the Cad protein gradient: (1) the gradient should be complementary to the anterior (Bcd) gradient; (2) as development progresses, there should be increasing repression in the abdominal region (∼50–80% A–P position); and (3) a posterior stripe should develop from time class 6 in C14A onwards at around 80% A–P position. Tll/Hkb: we replaced protein expression data with mRNA expression data for some of our reverse-engineering runs (see Results). Due to the relative constancy of tll and hkb expression patterns over time [43], [49], we created time-invariant expression profiles for these genes by averaging boundary positions across all cleavage cycles and time classes.

We use gene circuits for model fitting (reverse engineering) as described in [1], [38], [39], [46]–[51]. In brief, a gene circuit is a hybrid dynamical model incorporating discrete mitotic divisions of nuclei, as well as continuous gene regulatory dynamics within each nucleus. Each cleavage cycle consists of interphase, mitosis and division. During interphase, the change in gene product concentration (representing mRNA or protein, depending on the simulation) for gene a in nucleus i over time t is governed by the following of ordinary differential equation (ODE)(1)The three terms on the right-hand side of the equation define regulated gene product synthesis, diffusion and decay respectively. is the maximum synthesis rate; the diffusion rate (dependent on the distance between nuclei, which halves at every cleavage division: n defines the number of previous divisions); is the decay rate for the product of gene a. The sigmoid regulation-expression function captures the basic regulatory dynamics, and is defined as follows:(2)where(3)with the set of trunk gap genes defined as , and the set of external inputs as . Matrices W and E define the interactions between, respectively, the trunk gap genes themselves, and between the external inputs and trunk gap genes. The elements of these matrices, wba and ema, are called regulatory weights. These weights define the effect of b on a (or m to a) which can be (1) positive (activating gene product synthesis), (2) negative (inhibiting synthesis), or (3) (close to) zero (no regulatory interaction). ha is a threshold parameter representing uniformly distributed maternal factors. During mitosis, gene product synthesis is set to zero. After mitosis follows division, which is instantaneous. At division, gene product concentrations are copied equally to both daughter nuclei.

Finally, diffusion is implemented with no-flux boundary conditions.

Our models cover the trunk region of the embryo, from 35 to 87% A–P position. This region is somewhat reduced compared to 35 to 92% for earlier protein models [46]–[51], but covers the same set of gap domains, due to the slightly more anterior position of mRNA vs. protein domains. This results in gene circuits consisting of systems of 108 ODEs at C13, and 212 ODEs at C14A. Gap gene circuits were solved numerically from the beginning of C13 (t = 0 min) when gap proteins reach detectable levels, to the onset of gastrulation and the end of C14A (t = 71.100 min). We use the same division schedule as in [46], [48], [49]: mitosis occurs from 16.0 min to 21.0 min. At the end of mitosis, nuclear division takes place. Initial conditions for gap genes were calculated by interpolation between data points at C12 (t = −6.200 min) and C13 (t = 10.550 min) using the same temporal scaling scheme as described in the previous section. Initial conditions of the external inputs were taken from [49]. Time classes in C14A correspond to the following time points (in minutes): T1, t = 24.225; T2, t = 30.475; T3, t = 36.725; T4, t = 42.975; T5, t = 49.225; T6, t = 55.475; T7, t = 61.725; T8, t = 67.975 (see Figure 2 in [48]).

We follow a reverse engineering protocol as described in [46], [48], [49]. To estimate the values for parameters W, E, R, h, D and λ of the gene circuit model we performed global optimisation by means of parallel Lam Simulated Annealing (pLSA) [60]–[62] on the Mare Nostrum supercomputer at the Barcelona Computing Center (BSC; http://www.bsc.es). Per optimisation run we used 50 processor cores for an average duration of about 7 hours.

Simulated Annealing requires that candidate solutions have an associated cost (or energy) function that is minimised during the optimisation. We adapt the cost function from [49] as follows:(4)with T the set of time points (C13, C14A: T1–8) at which expression data is available, Nc(n) the number of nuclei after n divisions (50 at C13, 100 at C14A), representing positive weights associated with each data point, and referring to the expression level of gene product a at nucleus I and time t as derived from the experiments. If weights are all set to 1.0, the ‘cost’ equation represents a fit by ordinary least squares (OLS), which is the cost function previously used with gene circuit models [46]–[51], [63]. For fitting by weighted least squares (WLS), we use variable weights, which are inversely related to the level of expression (calculated as described in ‘Additional Data Processing for Model Fitting and Comparison’). This penalises ectopic gap gene expression, and improves frequency and quality of good fitting solutions as reported in [49].

Based on previous studies using gap gene circuit models, we fix certain parameters without negatively affecting the quality of the fits [46], [48], [49]. In the gene interaction matrix E we fix interactions of Hkb to zero, with the exception of Hkb→hb. Furthermore, we take ha = −2.5, for all gap genes , and for respectively [49].

To report the goodness of a fit, we use the root mean square (RMS), defined by(5)with the total number of data points in our data set. Since the RMS is independent of weights and the number of nuclei (data points) in the model, it allows us to compare WLS and OLS, as well as mRNA- and protein-based solutions. Solutions were selected for further analysis by several tests. Firstly, gene circuits were tested for numerical stability with respect to the solver (also known as solver sensitivity) and with respect to minute changes in parameter values (or the ‘brittleness’ of a solution). Subsequently, gene circuits were checked for visible gene expression patterning defects by means of visual inspection (see Text S1 for a categorisation of commonly encountered defects).