Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Robustness under Functional Constraint: The Genetic Network for Temporal Expression in Drosophila Neurogenesis

  • Akihiko Nakajima ,

    [email protected]

    Affiliation Department of Basic Science, University of Tokyo, Komaba, Tokyo, Japan

    ⨯
  • Takako Isshiki,

    Affiliation Center for Frontier Research, National Institute of Genetics, Mishima, Shizuoka, Japan

    ⨯
  • Kunihiko Kaneko,

    Affiliations Department of Basic Science, University of Tokyo, Komaba, Tokyo, Japan, ERATO Complex Systems Biology Project, Japan Science and Technology Agency, Komaba, Tokyo, Japan

    ⨯
  • Shuji Ishihara

    Affiliations Department of Basic Science, University of Tokyo, Komaba, Tokyo, Japan, PRESTO, Japan Science and Technology Agency, Kawaguchi, Saitama, Japan

    ⨯

Abstract

Precise temporal coordination of gene expression is crucial for many developmental processes. One central question in developmental biology is how such coordinated expression patterns are robustly controlled. During embryonic development of the Drosophila central nervous system, neural stem cells called neuroblasts express a group of genes in a definite order, which leads to the diversity of cell types. We produced all possible regulatory networks of these genes and examined their expression dynamics numerically. From the analysis, we identified requisite regulations and predicted an unknown factor to reproduce known expression profiles caused by loss-of-function or overexpression of the genes in vivo, as well as in the wild type. Following this, we evaluated the stability of the actual Drosophila network for sequential expression. This network shows the highest robustness against parameter variations and gene expression fluctuations among the possible networks that reproduce the expression profiles. We propose a regulatory module composed of three types of regulations that is responsible for precise sequential expression. This study suggests that the Drosophila network for sequential expression has evolved to generate the robust temporal expression for neuronal specification.

Author Summary

Cell fate specification is of key importance in the development of multicellular organisms. To specify various cell fates correctly, genetic networks precisely coordinate spatial and temporal gene expression patterns during various developmental stages. One central question in developmental biology is to elucidate the relationship between the pattern formation and the network architecture. During embryonic development of the Drosophila central nervous system, the neural stem cells express a group of genes in a definite order, which is responsible for the diversity of neural cells. To elucidate the underlying mechanism of the process, we analyzed the structure and dynamics of the genetic network for the temporal changes occurring in the Drosophila neural stem cells. Searching all the possible regulatory networks of these genes using a computer program, we detected the requisite regulations that reproduce observed gene expression profiles. By comparing the stability of the dynamics among the functional networks, we uncovered the robust nature of the actual Drosophila network against environmental and intrinsic fluctuations. These results indicate that the genetic network for sequential expression has evolved to be robust under functional constraints. Our study proposes regulatory modules that are responsible for the precise sequential expressions, which might exist in genetic networks for other temporal patterning processes.

Citation: Nakajima A, Isshiki T, Kaneko K, Ishihara S (2010) Robustness under Functional Constraint: The Genetic Network for Temporal Expression in Drosophila Neurogenesis. PLoS Comput Biol 6(4): e1000760. https://doi.org/10.1371/journal.pcbi.1000760

Editor: Nathan D. Price, University of Illinois at Urbana-Champaign, United States of America

Received: November 6, 2009; Accepted: March 24, 2010; Published: April 29, 2010

Copyright: © 2010 Nakajima et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Funding: This study was supported by a Grant-in-Aid for JSPS Fellows (21-10671 to AN) and by PRESTO of JST (to SI). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Precise coordination of cell fate decisions is crucial in the development of multicellular organisms. In the developmental processes, where a series of events occurs at a specific place and time, gene regulatory networks are responsible for implementing reliable biological functions [1], [2]. To obtain system-level understanding of such processes, it is necessary to integrate the molecular machinery of each regulation with architecture and dynamics at the regulatory network level. Biological functions achieved by gene networks are generally expected to possess robustness, i.e., insensitivity of system properties against a variety of perturbations that might originate from fluctuations during development and mutations through evolution. Recent investigations have addressed the questions of how robust biological functions are achieved through underlying molecular network architecture and its dynamic properties [3], [4], [5], [6], [7]. An illustrative example in developmental systems on this subject is segmentation of Drosophila melanogaster, which has been studied both experimentally and theoretically [8], [9], [10]. The requisite regulations or architecture of this system have been discussed at the network description level [10], [11], [12], [13], , and it is suggested that the underlying gene network has evolved to perform its processes in a robust manner [15], [16], [17].

Besides spatial patterning, temporal profiles of gene expression also play important roles in development [18], [19], [20]. Several computational studies have analyzed temporal expression profiles in biological processes such as the midgut development of sea urchin [21], [22] and vulval development of C. elegans [23]. These studies have shown relevant regulatory interactions and predicted unknown regulations for cell-fate specification.

The development of the Drosophila central nervous system (CNS) also manifests the importance of temporal patterning mechanism in development. Drosophila neural stem cell-like progenitors, called neuroblasts (NBs), generate a variety of neural cell types. During the embryonic development of the Drosophila CNS, NBs in the ventral nerve cord express certain transcription factors, i.e., Hunchback (Hb), Krüppel (Kr), Pdm1/Pdm2 (Pdm), and Castor (Cas), in a definite order (Fig. 1A–C) [24],[25],[26],[27]. In addition, the fifth factor, Seven-up (Svp), is expressed in the time window between Hb and Kr expression [28]. In association with this sequential expression, NBs divide asymmetrically to bud off a series of ganglion mother cells (GMCs). Each GMC undergoes an additional division to typically generate two postmitotic neurons. Depending on the transcription factors expressed in NBs at each division, postmitotic neurons acquire different cell fates. Thus, the sequentially expressed transcription factors control the cell-fate specification, thereby establishing the diversity of neurons in the Drosophila CNS. While neuronal specification process and generated cell types also depend on the spatial position [29], [30], [31] and lineage [32], [33] of NBs, the sequential expression is observed in a majority of ventral nerve cord NBs [34].

thumbnail
Download:
Figure 1. Sequential expression of temporal transcription factors within neuroblasts in the Drosophila CNS.

(A) The relative position of neuroblasts (NBs) in Drosophila embryo. The picture is the ventral view of NBs and shows Cas expression in the NBs at developmental stage 12. The bracket indicates a single segment. Dashed line corresponds to the midline. Scale bar: . (B) The expression levels of Hb, Kr, Pdm, and Cas in a single NB (NB 2–4 lineage) are shown from the developmental stage 10 to 12: early stage 10 (st. 10), early stage 11 (e11), mid stage 11 (11), late stage 11 (l11), mid stage 12 (12), late stage 12 (l12). (C) Schematic representation of the change of the expression pattern in a single NB. (D) The expression profiles of WT, loss-of-function, and overexpression mutants of the genes observed in the experiments (for references, see Table 2). (E) Reconstructed genetic network for sequential expression in Drosophila NBs. Repression from hb to cas (dashed line) was suggested to exist [26], although there is no direct verification. When the Drosophila network is invoked in this article, this regulation is also included. (F) Matrix representation of the Drosophila network.

https://doi.org/10.1371/journal.pcbi.1000760.g001

Isolated NBs exhibit sequential expression in vitro and differentiate into various neurons in a manner similar to that observed in vivo [35], [36]. Hb expression is switched off by Svp in a mitosis-dependent manner, while the subsequent expression of Kr, Pdm, and Cas proceeds in a mitosis-independent manner [28], [37]. These observations suggest that sequential expression of the genes is regulated cell-autonomously and occurs through mutual interactions among the factors.

In this study, we address the robustness of the gene network for sequential expression in the Drosophila CNS. One of the promising approaches to obtain insights into the system-level properties of biological systems is to compare the robustness of the actual network with that of other possible network architectures. Wagner considered how network architecture and robustness are related by studying circadian oscillation networks [38], although these networks lack a direct biological counterpart. Ma et al. studied the robustness of the Drosophila segmentation network [39], in which they had to arbitrarily eliminate components to reduce the size of the entire network. From theoretical and computational points of view, one advantage of studying temporal patterning in the Drosophila CNS is that the number of system components is so small that we can perform a comprehensive analysis of network architecture without any loss of biological relevance.

First, we explored the regulatory networks to reproduce the observed expression patterns in both wild-type (WT) and mutant embryos. We did not confine ourselves to only known regulations for sequential expression, but rather searched all possible networks that could reproduce the observed expression patterns. Studying the common structure of the specified genetic networks, we detected requisite regulations and predicted an unknown factor to reproduce known expression profiles. Second, we compared the robustness of the actual Drosophila network with that of the other networks reproducing the expression profiles. As a measure of robustness, we analyzed the stability of sequential expression against parameter variations and gene expression fluctuations. We found that the Drosophila network is highly robust and stable among possible functional networks. By further investigating the regulations necessary for the Drosophila network to be robust, we detected the responsible regulations. We propose a regulatory module composed of three kinds of regulations that is responsible for precise sequential expression of the Drosophila network.

Results

Temporal patterning network of D. melanogaster NBs

Expression profiles of temporal transcription factors (hb, Kr, pdm, cas, and svp) in Drosophila NBs are summarized in Figure 1D for WT, loss-of-function, and overexpression embryos [25], [26], [28], [36], [40], [41]. It has been considered that these sequential expressions are produced (or at least modulated) by mutual regulations among the temporal transcription factors [24], [25]. We reconstructed the gene network for sequential expression in Drosophila NBs from the literature as shown in Figure 1E and F (for references, see Table 1).

thumbnail
Download:
Table 1. List of the regulatory interactions of the genes in the NB temporal patterning network.

https://doi.org/10.1371/journal.pcbi.1000760.t001

Modeling gene network dynamics by Boolean description

First, we searched for regulatory networks that reproduce the sequential expression patterns of both WT and mutants. To investigate gene expression dynamics, we adopted a Boolean-type model [6] (see Materials and Methods for details of the model and the following analysis):(1)where represents the expression state of gene i () at the t-th time step and takes either 1 (ON) or 0 (OFF). Regulation from gene j to gene i is either positive (Jij >0), negative (Jij <0), or zero (Jij = 0), which corresponds to activation, repression, or absence of regulation, respectively. The state of gene i at the next step () is 1 when the sum of regulatory inputs is positive () or 0 when the sum is negative (). When the sum equals zero (), takes the default expression state : . In this study, the value of Jij is supposed to take one of the discrete values . The large negative value (−5) of Jij signifies that the expression of a gene is completely shut off in the presence of a repressor. This choice of large negative value comes from experimental observations of mutants. In experimentally observed expression patterns (Fig. 1D), genes are not activated when both repressors and activators are expressed. For example, in Kr++ and pdm++ embryo (here “++” means overexpression of the gene), pdm and cas expression is not observed in hb-expressing time window, although their activators are overexpressed. This indicates that the repressive effect from hb is dominant over pdm activation by Kr and cas activation by pdm.

Initial expression state of genes is set to 0, except for Hb, which emulates the NB gene expression in the first stage of sequential expression [24], [25]. Thus far, the only known function of Svp during the early stage is downregulation of Hb. There is no evidence that Svp regulates or is regulated by other temporal transcription factors during the expression series: Kr Pdm Cas [28]. In addition, Hb is only regulated by Svp and not by the other three factors (Kr, Pdm, and Cas). Thus, in the model, we assumed a pulsed expression of Svp as an input to the system, resulting in downregulation of Hb at the next time step. The temporal expression dynamics of Kr, Pdm, and Cas follow Eq. (1) with assigned values of Jij (Fig. 1F).

The regulatory networks of known factors do not reproduce the experiments

Based on the above formulation, we investigated whether the reconstructed Drosophila gene network (Fig. 1E and F) is sufficient to reproduce the sequential expression observed in WT, as well as all the known single loss-of-function and overexpression mutants, i.e., hb−, Kr−, pdm−, cas−, hb++, Kr++, pdm++, and cas++ (Fig. 1D, Table 2). Presently, we cannot specify the value of the parameters , and from empirical data; thus, each value could be arbitrarily chosen from (). We studied all 23 combinations of and found that the dynamics coincide with the expression profile in WT but not in some mutants for each choice of parameters (examples shown in Fig. 2). Depending on the parameter values, the expression dynamics changed to some extent, but none of the possible combinations reproduced the expression profiles of all of the mutants. For example, in case of , , and , the dynamics of the network for hb− and Kr− did not agree with the experiments (Fig. 2A), and in case of , , and , the dynamics of hb− and pdm− did not (Fig. 2C).

thumbnail
Download:
Figure 2. Reconstructed Drosophila network cannot reproduce the experimentally reported expression profiles.

Sequential gene expression of reconstructed Drosophila network is simulated using Boolean model. The grids filled with colors represent ON states of the genes. The dynamics could be different depending on the choice of the default expression states . (A) ; (B) ; and (C) .

https://doi.org/10.1371/journal.pcbi.1000760.g002

thumbnail
Download:
Table 2. List of references for the sequential expression pattern in various genotypes.

https://doi.org/10.1371/journal.pcbi.1000760.t002

We then investigated whether networks other than the Drosophila network can reproduce the observed expression profiles by checking all the 312 ( = 531,441) combinations of Jij values. The dynamics agreed with the expression profile in WT for a large number of networks (39,391 out of 531,441), but any networks composed of hb, Kr, pdm, cas, and svp did not reproduce the profiles in both WT and mutants.

Introduction of a presumptive factor is sufficient to reproduce the expression profiles

Preceding results indicate the difficulty of reproducing the observed expression patterns only with known constituents. We therefore introduced an additional presumptive regulator (x). The expression state of x was assumed to start in the ON state and change into OFF, or vice versa at () (see Materials and Methods). Including this assumption, we reinvestigated the dynamics of all 315 ( = 14,348,907) possible regulatory networks with all the possible switching timings of x. In the case that the expression of x switches OFF to ON, none of the networks conformed to the expected expression profiles. On the other hand, in the case that the expression of x switches ON to OFF, we found that 384 networks (<0.003%) reproduced the expression profiles of both WT and mutants. We refer to the detected networks as “the functional networks” hereafter in the study.

Comparing the regulatory interactions of the functional networks, we found that the regulations shared among all the functional networks are coincident with experimentally verified regulations (colored as black in Fig. 3A). In addition, activation of Kr and repression of cas by a presumptive factor x appear in all of the functional networks (colored as brown in Fig. 3A). The genetic network composed of these common regulations is a minimum network to reproduce the expression profiles of WT and mutants. To quantify the similarity among the functional networks, we measured the distances of the 384 functional networks from the actual Drosophila network (Fig. 3C); the distances are defined by the number of different regulations (see Materials and Methods). As a reference, we also performed the same analyses of distance measurement for all possible networks and the networks that are randomly reconnected from functional networks (see Materials and Methods). For all possible networks, the frequency distribution of the distances shows that the network architectures are different from the actual Drosophila network by 7.81.5 regulations. The reconnected networks yield similar results, albeit with slightly decreased distances (7.01.7 regulations). In contrast, the architectures of the functional networks differ by only 2.41.1 regulations. The architectures of the functional networks resemble that of the actual Drosophila network. These indicate that the gene networks that reproduce the known sequential expression patterns are highly constrained in their topologies.

thumbnail
Download:
Figure 3. Architecture of the detected functional networks.

(A) Architecture of the functional networks reproducing the gene expression profiles observed in the experiments. The black arrows are the regulations that appear in all the functional networks. The brown arrows are the regulations from the presumptive factor x that also appear in all the functional networks. The other regulations existing in the actual Drosophila network are shown by gray arrows. (B) Matrix representation of the functional networks. Elements of {Jij} are shown as either + for activation, − for repression, or 0 for the absence of regulation. (C) Frequency distributions of the distances of networks from the Drosophila network. The distributions are drawn from the functional networks (N = 384; magenta), all the possible networks (N = 14,348,907; blue), and the networks randomly reconnected from the functional ones (N = 38,400; yellow). From each of the functional networks, 100 reconnected networks were generated. The regulatory interactions from x and positive self-feedbacks are neglected in counting the number of different regulations.

https://doi.org/10.1371/journal.pcbi.1000760.g003

Robustness of the Drosophila network against parameter variations and expression noise

Because there are multiple network architectures that explain the observed expression profiles as shown above, we then investigated the characteristics of the actual Drosophila network among the functional networks. From the biological point of view, the sequential expression in NBs should proceed reliably despite developmental disturbances such as cell-to-cell variation and intracellular fluctuations. We thus evaluated the stability of sequential expression for each of the detected functional networks and compared the properties of the actual Drosophila network to those of the other networks. To address the problem quantitatively, we extended the previous Boolean model into a model of ordinary differential equations with fluctuations in gene expression, where the concentrations of mRNAs {Mi(t)} and proteins {Pi(t)} obey the following equations [42], [43] (see Materials and Methods for the details of the model and the following analysis):(2)Here i refers to one of each gene: . The variables {Mi(t)} and {Pi(t)} take continuous values, unlike the previous Boolean description. The precise function form of promoter activities {Fi({Pj(t)})} is dependent on the regulatory interactions of the genetic networks and the default promoter activities {Si}, corresponding to the Boolean model. The time-dependent variables represent the noise in promoter activities. Here we have assumed that the expression noise comes from the transcription process (noise is incorporated only in the dynamics of {Mi(t)}). One reason is the practical convenience in the numerical calculations. In addition, recent quantitative analyses of gene expression have indicated that the gene expression noise mainly arises from transcription [44], [45], [46]. However, we should note that the result and conclusion obtained from the following analysis does not change even if we incorporate noise in the dynamics of {Pi(t)} as well (data is not shown).

Typical dynamics of the Drosophila network are shown in Figure 4, where sequential expression of WT is reproduced. The dynamics of the model are largely dependent on the parameter values and the noise intensities, and coincide with the experimental observations only under appropriate conditions. Therefore, such sensitivity to parameter variation is important for the development to proceed under environmental and individual fluctuations.

thumbnail
Download:
Figure 4. Temporal dynamics of the Drosophila network in the continuous model.

The dynamics of expression levels of proteins {Pi(t)} with different parameter values (upper) and discretized representation of a typical temporal dynamics (lower). In addition to the known genes, the presumptive factor x is also incorporated. The expression level of X changes from a high level to a low level as in the previous model. Each gene is considered to be in the ON state when the expression level is larger than a threshold Pth. The parameter values of and {Si} are randomly selected from the following ranges: for and for ; and and . The other parameter values are set as shown in Table 4.

https://doi.org/10.1371/journal.pcbi.1000760.g004

To characterize sensitivity, we measured the fraction of successes; that is, the fraction of the parameter sets that can reproduce the expression profile of WT among all the trials of random parameter assignments [15], [39]. To judge whether the dynamics coincide with the expression profile in Drosophila NBs, the dynamics of the protein concentrations {Pi} were discretized to 1 (0) for Pi > Pth (Pi < Pth). The threshold Pth was set as Pth = 0.2. The temporal dynamics of a network were accepted when the discretized dynamics satisfied the condition for WT in Table 3. To obtain the effect of parameter variation, we carried out the simulation without stochastic terms in Eq. (2). In each network, we repeated the simulations with random assignment of parameter values and calculated the fraction of successes (Fig. 5A). The Drosophila network scored the highest fraction of successes among the functional networks, and the networks closer to the Drosophila network tended to have higher scores.

thumbnail
Download:
Figure 5. Robustness of the gene expression profiles in the functional networks.

(A) The fraction of trials that reproduce the experimental expression profile against random assignments of parameters. The values of , , and are randomly chosen within the ranges shown in Table 4. The other fixed parameter values are also listed in Table 4. Neglecting the positive self-feedback regulations in the 384 functional networks, 120 networks were chosen and investigated (Materials and Methods). The dynamics were checked for 50,000 trials in each network. The networks were sorted based on the distance from the Drosophila network (Nd). Here Nd corresponds to the number of regulations different from the Drosophila network. Because there are a few possible regulations from the unknown factor x, more than one network with Nd = 0 exist. (B) The fractions of the trials that reproduce the experimental profile under expression noise (vertical axis) are plotted against the fraction of successes against the random parameter assignments. To analyze the stability against noise, we used 1000 different parameter sets, by which the expression profile is reproduced in the absence of noise for each network. The dynamics were checked for 50 trials for each parameter set.

https://doi.org/10.1371/journal.pcbi.1000760.g005

thumbnail
Download:
Table 3. Criterion for expression profile in each genotype.

https://doi.org/10.1371/journal.pcbi.1000760.t003

We also investigated the dynamical stability of the gene networks against fluctuations. In this case, we performed the stochastic simulations in Eq. (2) with expression noise. To evaluate stability against noise, we chose the parameter values with which the expression profile is reproduced in the absence of noise. We then measured the relative fraction of successes under fluctuation. As is shown in Figure 5B, the fraction of successes under expression noise increased with the similarity to the actual Drosophila network as the fraction of successes under parameter variations. Thus, the Drosophila network lies at the top level of the functional networks in terms of robustness against these perturbations.

Regulations that heighten functional stability

Because the Drosophila network has several other regulations in addition to the minimum functional network (gray arrows in Fig. 3A), these regulations might be responsible for the robustness shown above. We compared the robustness among the networks with or without the additional regulations. The fraction of successes against parameter variations for these networks is plotted in Figure 6A. The minimum network reproduces the sequential expression under the appropriate parameters, but the robustness is much lower than that of the Drosophila network. The scores of networks that lack one of the regulations fall between the minimum and the Drosophila network. Stability to expression noise was also evaluated by changing noise intensity, and similar results were obtained (Fig. 6B). The fraction of successes decreased as the noise intensity became larger, but the effect of noise on the Drosophila network was less severe than that on the minimum network. Thus, each of these regulations contributes to the robustness of the system.

thumbnail
Download:
Figure 6. Contribution of the actual regulations to the robustness of the system.

(A) The fraction of the trials that reproduce the experimental WT expression against parameter variations. The data of Figure 5A are replotted for the Drosophila network, the networks lacking an indicated regulation (one of the gray arrows in Fig. 3A) and the minimum network (black and brown arrows in Fig. 3A). (B) The fractions of the trials that reproduce the experimental profile under gene expression noise with various intensities. We used 5,000 different parameter sets with which the profile is reproduced in the absence of noise. The dynamics are checked for 50 trials for each parameter set.

https://doi.org/10.1371/journal.pcbi.1000760.g006

To elucidate the roles of these regulations, we tried random parameter assignments for each of these networks and sampled successful parameter sets that reproduce WT sequential expression profile (Fig. 7). In the Drosophila network (Fig. 7A), wide ranges of parameter values are allowed, indicating that this network reproduces the required profile without quantitative tuning of parameters, and thus, shows high robustness. For other networks (Fig. 7B–E), the ranges are narrower for some parameters (as clearly seen in Spdm and Scas), and the numbers of successful parameter sets are less than those obtained for the Drosophila network.

thumbnail
Download:
Figure 7. Graphical representation of parameter sets with which the WT sequential expression profile is reproduced.

(A) The Drosophila network, the networks lacking (B) activation from Kr to pdm, (C) activation from pdm to cas, (D) repression from hb to cas, and (E) the minimum network. The parameters involved in minimum network are shown. Each spoke represents a value of indicated parameter between the range used for random parameter assignment (Table 4). The value of is shown by normal scale and those of the other parameters are shown by log scale. Each polygon indicates one parameter set. Solid and broken lines indicate mean and s.d. of obtained parameters. The data are drawn from 5,000 trials of the random assignment of parameter values.

https://doi.org/10.1371/journal.pcbi.1000760.g007

How is the robust nature of the Drosophila network implemented by these regulations? As seen above, the parameter values of Spdm and Scas (default promoter activities of pdm and cas) are most influenced by the loss of these regulations. Because expression of a gene is induced by either the activity of the default promoter or the activators (see Materials and Methods), additional regulations in the Drosophila network (gray arrows in Fig. 3A) might compensate for the loss of default activities. To verify this possibility, we measured the dependence of the fraction of successes on the strength of regulations (, , and ) and default promoter activities (Spdm and Scas) (Fig. 8A–C).

thumbnail
Download:
Figure 8. Parameter dependencies of robustness for the Drosophila network.

The fractions of successes for random assignment of parameter values are plotted under the different strengths of regulations (, , and ) and default promoter activities (Spdm and Scas). Dependencies of robustness to (A) (strength of activation from Kr to pdm) and Spdm, (B) (strength of activation from pdm to cas) and Scas, and (C) (strength of the repression from hb to cas) and Scas. The other parameters are set as listed in Table 4. The temporal dynamics were tested for 50,000 trials.

https://doi.org/10.1371/journal.pcbi.1000760.g008

Figure 8A shows the fraction of successes for random assignments of parameter values under given strengths of and Spdm. To score high reproducibility, Spdm must be large for small , but need not to be large for sufficiently large . This indicates that activation of pdm expression by Kr indeed compensates for the loss of default promoter activity of pdm. Thus, for the network lacking this regulation, the default promoter activity is necessary because inductions from other factors are absent. A similar relationship is found between and Scas (Fig. 8B).

As for repression of cas by hb, the role for robustness seems to be different from the above two. When the absolute value of is small, Scas must be small to achieve a high fraction of successes (Fig. 8C). As becomes larger, a higher value of Scas is allowed. This is because the repression from hb to cas reduces the mis-expression of cas in the early stage of sequential expression. Grosskortenhaus et al. suggested the direct repression from hb to cas [26], although there is no confirmative evidence to our knowledge. This regulation possibly contributes to the robustness of the actual system.

Discussion

Through the present analyses, we obtained 384 functional networks that reproduce the sequential expression of both WT and mutants. The detected functional networks exhibit high similarity in regulatory interactions among the transcription factors (Fig. 3). This exemplifies the importance of the regulations in the minimum network for the sequential expression. In addition, the actual Drosophila network scores quite high on reproducibility of the WT sequential expression among all the functional networks (Fig. 5 and 6). Below, we discuss the biological implications of the temporal patterning of Drosophila NBs drawn from our numerical analyses.

Existence of an unknown factor can reproduce the expression patterns of WT and mutants

In this study, we introduced an additional presumptive factor x to obtain networks that reproduce the sequential expression of both WT and mutants. Because x is hypothetical, we discuss its validity here.

Because the loss-of-function mutant of any one gene has only minor effects on the expression sequence (Fig. 1D), several previous reports suggested the existence of either unknown regulators or an additional clock mechanism that regulates the sequential expression [25], [26]. Our assumption is feasible for explaining experimental results because it does not need any other clock mechanism or superfluous multiple regulators. It is notable that our analysis indicates that the possible regulations of the presumptive factor are highly restricted; the expression of x switches ON state to OFF state (Fig. 4), and all the functional networks have activation of Kr and repression of cas by x (Fig. 3A). Thus, our assumption can be tested in future experiments in vivo.

We should note that while the regulator x is needed to explain the mutant profiles under our modeling assumptions, the mutual regulations of only known factors also reproduce the WT sequential expression (Fig. 1D). Therefore, the regulations among hb, Kr, pdm, and cas would play a primary role as discussed below.

Minimum network for the sequential expression

An effective way to capture network function is to focus on the specific substructures (network motifs or modules) [1], [13], [14], [16], [39], [47]. Comparing all the functional networks, we detected the minimum structure for the sequential expression, which contains two successive regulatory loops (Fig. 3A and 9A); one is composed of hb, Kr, and pdm, and the other of Kr, pdm, and cas. In each loop, one gene represses the previous and the second next factor. The repressions of the second next factors (hb to pdm and Kr to cas) define the induction timing of the regulated factors, since they are kept repressed until the regulators are switched off. The feedback repression of the previous factors (pdm to Kr and cas to pdm) ensures their downregulation, which promotes the progress of the sequential expression. These coincide with the observations by Kambadur et al., who experimentally showed that the repressions from hb and cas define the temporal window of Pdm [24]. These repressive regulations and the activation from hb to Kr compose the minimum network for sequential expression (Fig. 9A). Although they are enough to reproduce the sequential expression under appropriate conditions, the expression profiles could be easily perturbed by parameter variations or increase of noise (Fig. 5 and 9A).

thumbnail
Download:
Figure 9. Regulatory module for precise sequential expression.

The regulatory interactions and schematic expression profiles of the networks. (A) Regulatory interactions of the minimum network for sequential expression (left). This network reproduces the sequential expression under appropriate conditions (middle). However, the parameter variations from the appropriate values and the increase of noise could easily alter the expression profiles (right). (B) Regulatory interactions of the Drosophila network (left). Three types of regulations in this network enable the temporal expression in the precise order.

https://doi.org/10.1371/journal.pcbi.1000760.g009

Robustness of the Drosophila network: mechanism generating the precise sequential expression

In the two loops of the Drosophila network, activations from one gene to the next (Kr to pdm and pdm to cas) exist in addition to the repressive regulations. Other functional networks do not necessarily have these activations, but the activations can compensate for the loss of default promoter activities (Fig. 8A and B). These regulations achieve precise expression by enhancing the correlations among the factors and heightening the stability against fluctuations (Figs. 5B and 6B). From these results, we conclude that three types of regulations (activation of the next factor, feedback repression, and repression of the second next factor) compose a regulatory module for precise temporal expression, as summarized in Figure 9B. The feature of this network module embodies the robustness of the Drosophila network.

Do the previous discussions have any implications on other developmental processes? In the studies of spatial patterning in Drosophila segmentation, it was claimed that the frequent substructure feed forward loop (FFL) can set the positions of expression domains [13], and mutual feedback repressions between the gap genes also have a pivotal role in the formation of expression domains with steep boundaries [12], [47]. In case of the Drosophila network for sequential expression, preceding genes activate the next ones, while these genes repress the preceding ones. Similar regulatory interactions are reported in the yeast cell cycle by Lau et al. [48]. Thus, such asymmetric mutual regulations would be a general mechanism that serves as precise switches in the process of temporal patterning.

Role of the robustness in Drosophila neurogenesis

We showed that the temporal specification network of Drosophila NBs contains not only the regulations necessary for generating sequential expression, but also additional regulations to achieve higher precision in the expression. In each hemisegment of Drosophila embryo, 30 different NBs are generated through spatial heterogeneity [29]. To guarantee sequential expression of common temporal transcription factors despite their differences in Drosophila NBs, the robustness of the system might be important.

The robust nature of the Drosophila temporal network could be the consequence of evolutionary optimization in the reproducibility of the sequential expression under functional constraint. In future, we expect that experimental manipulation of corresponding enhancers will be able to clarify the relevance of each regulation to temporal patterning and stability.

Materials and Methods

Analysis of temporal dynamics of the genetic networks with the Boolean model

Here we describe the details of the Boolean model (Eq. (1)). The expressions of svp and x occur as inputs to the system. A pulse of svp expression always occurs at t = 1. Expression of x switches either from ON to OFF state, or from OFF to ON state at (). Once we assign the switching time of x expression (), its value becomes fixed through the analysis of expression patterns for all the genotypes. Because the autonomous pulsed expression of svp results in hb downregulation, we set Jhb, svp = −5, Jhb, j = 0 (j = hb, Kr, pdm, cas, or x), and Jk, svp = 0 (k = Kr, pdm, or cas) throughout this study. The time step at which we finish the simulation () was set as .

We thus investigated the behaviors of the remaining three factors (Kr, pdm, and cas) under the given regulatory interactions {Jij}. The total number of combinations of the parameters is 3M23 (the number of possible network architecture {Jij} multiplied by the number of default expression states ), where M is the number of regulations. To simulate the dynamics for mutants, we always set the expression state of the corresponding gene to 0 (OFF) for loss-of-function or to 1 (ON) for overexpression. We then examined whether the temporal dynamics of the genetic networks are coincident with the expression profiles of each mutant (Fig. 1D and Table 3).

Analysis of network statistics

In order to measure the similarity between the functional networks and the actual Drosophila network, we used two types of network ensembles as references. One is the ensemble of the possible network architectures. The other is a set of reconnected networks generated from the functional networks by iterative random reconnections of the matrix elements (1,000 iterations). The numbers of positive and negative regulations are preserved in the iterations.

To count the number of different regulations between functional networks and the actual Drosophila network, we neglected the regulations from x and positive self-feedbacks because the existence of those is uncertain from the experimental data.

Continuous model of the expression dynamics

We introduced the continuous model with stochasticity as shown in Equation (2). The promoter activity of gene i (i = hb, Kr, pdm, cas, or x) is described as follows,Regulatory interactions are continuous equivalents of {Jij} in the Boolean model, and g(x) is a piece-wise linear function such that g(x) = x for x>0 and g(x) = 0 for x<0. The parameters {Si} are the default activities of the promoters. Transcription of a gene is induced when the total regulatory inputs become positive (), and is suppressed when they become negative (). In order to consider the effect of fluctuations on the expression dynamics, we introduced additive white Gaussian noise : (Eq (2)), where is the noise intensity of gene i.

The expression of hb and x is induced only by the default promoter activities because all the regulations are absent for these two (). To describe the expression change of hb and x, the promoter activities of these two are set as Shb >0 for (Sx >0 for ) and Shb = 0 for (Sx = 0 for ), respectively. The promoter activities of the others are always assumed to exist (SKr, Spdm, and Scas >0). The noise intensities are also set as (>0) for and for (i = hb, x). Those of the other genes are (>0) (j = Kr, pdm, cas), Here we simply assume that the noise intensities of the genes take the same value . The noise intensity is set as in Figure 4, and in Figure 5. Noise intensity (horizontal axis) in Figure 6B means the value of .

Analysis of the robustness of the networks

For the continuous model, we considered two different types of robustness: (1) the reproducibility of the sequential expression against parameter variations and (2) dynamical stability against temporal fluctuations. To analyze the former, the default promoter activities {Si} were assigned randomly within the defined ranges. The values of the matrix were set to 0 when the corresponding regulations were absent (the corresponding element of the Boolean model takes Jij = 0) or assigned randomly when they are present (Jij0). In order to confine our attention to the properties of network architectures, the other parameters (, , , and ) were fixed throughout the analysis. The ranges and the fixed values of the parameters are listed in Table 4. Robustness against temporal fluctuations is measured as explained in the main text.

thumbnail
Download:
Table 4. Parameter values used for continuous dynamics of the genetic networks.

https://doi.org/10.1371/journal.pcbi.1000760.t004

In the simulations, we found that the existence of positive self-regulation enhanced the fraction of successes in many cases, but hardly affected the sequential expression. To focus on the contributions of mutual regulations of genes to robustness, we neglected the positive self-feedback regulations and confined the analysis to 120 out of 384 functional networks.

Acknowledgments

We thank T. Suzuki and K. Fujimoto for their useful comments on the manuscript, and H. Takagi and T. Shibata for their discussion.

Author Contributions

Conceived and designed the experiments: AN SI. Performed the experiments: AN. Analyzed the data: AN. Wrote the paper: AN TI KK SI.

References

  1. 1. Alon U (2006) An introduction to systems biology: design principles of biological circuits: Chapman & Hall.
  2. 2. Levin M, Davidson E (2005) Gene regulatory networks for development. Proc Natl Acad Sci USA 102: 4936–4942.
  3. 3. Wagner A (2005) Robustness and evolvability in living systems: Princeton University Press.
  4. 4. Ciliberti S, Martin O, Wagner A (2007) Robustness can evolve gradually in complex regulatory gene networks with varying topology. PLoS Comput Biol 3: e15.
  5. 5. Kaneko K (2007) Evolution of robustness to noise and mutation in gene expression dynamics. PLoS One 2: e434.
  6. 6. Li F, Long T, Lu Y, Ouyang Q, Tang C (2004) The yeast cell-cycle network is robustly designed. Proc Natl Acad Sci USA 101: 4781–4786.
  7. 7. Barkai N, Leibler S (1997) Robustness in simple biochemical networks. Nature 387: 913–917.
  8. 8. Akam M (1987) The molecular basis for metameric pattern in the Drosophila embryo. Development 101: 1–22.
  9. 9. Schroeder M, Pearce M, Fak J, Fan H, Unnerstall U, et al. (2004) Transcriptional control in the segmentation gene network of Drosophila. PLoS Biol 2: E271.
  10. 10. Perkins T, Jaeger J, Reinitz J, Glass L (2006) Reverse engineering the gap gene network of Drosophila melanogaster. PLoS Comput Biol 2: e51.
  11. 11. Sánchez L, Thieffry D (2003) Segmenting the fly embryo: a logical analysis of the pair-rule cross-regulatory module. J Theor Biol 224: 517–537.
  12. 12. Jaeger J, Blagov M, Kosman D, Kozlov K, Manu , et al. (2004) Dynamical analysis of regulatory interactions in the gap gene system of Drosophila melanogaster. Genetics 167: 1721–1737.
  13. 13. Ishihara S, Fujimoto K, Shibata T (2005) Cross talking of network motifs in gene regulation that generates temporal pulses and spatial stripes. Genes Cells 10: 1025–1038.
  14. 14. Fujimoto K, Ishihara S, Kaneko K (2008) Network evolution of body plans. PLoS One 3: e2772.
  15. 15. von Dassow G, Meir E, Munro E, Odell G (2000) The segment polarity network is a robust developmental module. Nature 406: 188–192.
  16. 16. Ingolia N (2004) Topology and robustness in the Drosophila segment polarity network. PLoS Biol 2: e123.
  17. 17. Manu , Surkova S, Spirov A, Gursky V, Janssens H, et al. (2009) Canalization of gene expression and domain shifts in the Drosophila blastoderm by dynamical attractors. PLoS Comput Biol 5: e1000303.
  18. 18. Thummel C (2001) Molecular mechanisms of developmental timing in C. elegans and Drosophila. Dev Cell 1: 453–465.
  19. 19. Pearson B, Doe C (2004) Specification of temporal identity in the developing nervous system. Annu Rev Cell Dev Biol 20: 619–647.
  20. 20. Caygill E, Johnston L (2008) Temporal regulation of metamorphic processes in Drosophila by the let-7 and miR-125 heterochronic microRNAs. Curr Biol 18: 943–950.
  21. 21. Yeo Z, Wong S, Arjunan S, Piras V, Tomita M, et al. (2007) Sequential logic model deciphers dynamic transcriptional control of gene expressions. PLoS One 2:
  22. 22. Yuh C, Bolouri H, Davidson E (2001) Cis-regulatory logic in the endo16 gene: switching from a specification to a differentiation mode of control. Development 128: 617–629.
  23. 23. Fisher J, Piterman N, Hajnal A, Henzinger T (2007) Predictive modeling of signaling crosstalk during C. elegans vulval development. PLoS Computational Biology 3: e92.
  24. 24. Kambadur R, Koizumi K, Stivers C, Nagle J, Poole S, et al. (1998) Regulation of POU genes by castor and hunchback establishes layered compartments in the Drosophila CNS. Genes Dev 12: 246–260.
  25. 25. Isshiki T, Pearson B, Holbrook S, Doe C (2001) Drosophila neuroblasts sequentially express transcription factors which specify the temporal identity of their neuronal progeny. Cell 106: 511–521.
  26. 26. Grosskortenhaus R, Robinson K, Doe C (2006) Pdm and Castor specify late-born motor neuron identity in the NB7-1 lineage. Genes Dev 20: 2618–2627.
  27. 27. Novotny T, Eiselt R, Urban J (2002) Hunchback is required for the specification of the early sublineage of neuroblast 7-3 in the Drosophila central nervous system. Development 129: 1027–1036.
  28. 28. Kanai M, Okabe M, Hiromi Y (2005) seven-up Controls switching of transcription factors that specify temporal identities of Drosophila neuroblasts. Dev Cell 8: 203–213.
  29. 29. Technau G, Berger C, Urbach R (2006) Generation of cell diversity and segmental pattern in the embryonic central nervous system of Drosophila. Dev Dyn 235: 861–869.
  30. 30. Egger B, Chell J, Brand A (2008) Insights into neural stem cell biology from flies. Philos Trans R Soc Lond B Biol Sci 363: 39–56.
  31. 31. Wheeler S, Stagg S, Crews S (2008) Multiple Notch signaling events control Drosophila CNS midline neurogenesis, gliogenesis and neuronal identity. Development 135: 3071–3079.
  32. 32. Baumgardt M, Miguel-Aliaga I, Karlsson D, Ekman H, Thor S (2007) Specification of neuronal identities by feedforward combinatorial coding. PLoS Biol 5: e37.
  33. 33. Udolph G, Rath P, Tio M, Toh J, Fang W, et al. (2009) On the roles of Notch, Delta, kuzbanian, and inscuteable during the development of Drosophila embryonic neuroblast lineages. Dev Biol 336: 156–168.
  34. 34. Karcavich R (2005) Generating neuronal diversity in the Drosophila central nervous system: a view from the ganglion mother cells. Dev Dyn 232: 609–616.
  35. 35. Brody T, Odenwald W (2000) Programmed transformations in neuroblast gene expression during Drosophila CNS lineage development. Dev Biol 226: 34–44.
  36. 36. Grosskortenhaus R, Pearson B, Marusich A, Doe C (2005) Regulation of temporal identity transitions in Drosophila neuroblasts. Dev Cell 8: 193–202.
  37. 37. Mettler U, Vogler G, Urban J (2006) Timing of identity: spatiotemporal regulation of hunchback in neuroblast lineages of Drosophila by Seven-up and Prospero. Development 133: 429–437.
  38. 38. Wagner A (2005) Circuit topology and the evolution of robustness in two-gene circadian oscillators. Proc Natl Acad Sci USA 102: 11775–11780.
  39. 39. Ma W, Lai L, Ouyang Q, Tang C (2006) Robustness and modular design of the Drosophila segment polarity network. Mol Syst Biol 2: 70.
  40. 40. Cleary M, Doe C (2006) Regulation of neuroblast competence: multiple temporal identity factors specify distinct neuronal fates within a single early competence window. Genes Dev 20: 429–434.
  41. 41. Tran K, Doe C (2008) Pdm and Castor close successive temporal identity windows in the NB3-1 lineage. Development 135: 3491–3499.
  42. 42. McAdams H, Arkin A (1998) Simulation of prokaryotic genetic circuits. Annu Rev Biophys Biomol Struct 27: 199–224.
  43. 43. Smolen P, Baxter D, Byrne J (2000) Modeling transcriptional control in gene networks–methods, recent results, and future directions. Bull Math Biol 62: 247–292.
  44. 44. Blake W, Kaern M, Cantor C, Collins J (2003) Noise in eukaryotic gene expression. Nature 422: 633–637.
  45. 45. Raser J, O'Shea E (2005) Noise in gene expression: origins, consequences, and control. Science 309: 2010–2013.
  46. 46. Newman J, Ghaemmaghami S, Ihmels J, Breslow D, Noble M, et al. (2006) Single-cell proteomic analysis of S. cerevisiae reveals the architecture of biological noise. Nature 441: 840–846.
  47. 47. Ishihara S, Shibata T (2008) Mutual interaction in network motifs robustly sharpens gene expression in developmental processes. J Theor Biol 252: 131–144.
  48. 48. Lau K, Ganguli S, Tang C (2007) Function constrains network architecture and dynamics: a case study on the yeast cell cycle Boolean network. Phys Rev E 75: 051907.