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

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.


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,14], 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].
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 mitosisdependent 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).

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): where X t i represents the expression state of gene i (i [ hb, Kr, pdm, cas f g ) at the t-th time step and takes either 1 (ON) or 0 (OFF). Regulation from gene j to gene i is either positive (J ij .0), negative (J ij ,0), or zero (J ij = 0), which corresponds to activation, repression, or absence of regulation, respectively. The state of gene i at the next step (X tz1 i ) is 1 when the sum of regulatory inputs is positive ( P j J ij X t j w0) or 0 when the sum is negative ( P j J ij X t j v0). When the sum equals zero ( P j J ij X t j~0 ), X tz1 i takes the default expression state h i : h i [ f0, 1g. In this study, the value of J ij is supposed to take one of the discrete values J ij [ f1, 0, {5g. The large negative value (25) of J ij 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,

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.  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.

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 2 , Kr 2 , pdm 2 , cas 2 , hb ++ , Kr ++ , pdm ++ , and cas ++ (Fig. 1D, Table 2). Presently, we cannot specify the value of the parameters h Kr , h pdm and h cas from empirical data; thus, each value could be arbitrarily chosen from h i [ f0,1g (i [ fKr, pdm, casg). We studied all 2 3 combinations of fh i g 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 h Kr~0 , h pdm~0 , and h cas~1 , the dynamics of the network for hb 2 and Kr 2 did not agree with the experiments ( Fig. 2A), and in case of h Kr~1 , h pdm~1 , and h cas~1 , the dynamics of hb 2 and pdm 2 did not (Fig. 2C).
We then investigated whether networks other than the Drosophila network can reproduce the observed expression profiles by checking all the 3 12 ( = 531,441) combinations of J ij 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 t~t switch (0ƒt switch ƒt end) ( see Materials and Methods). Including this assumption, we reinvestigated the dynamics of all 3 15 ( = 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.8+1.5 regulations. The reconnected networks yield similar results, albeit with slightly decreased distances (7.0+1.7 regulations). In contrast, the architectures of the functional networks differ by only 2.4+1.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.

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 quantita- Table 1. List of the regulatory interactions of the genes in the NB temporal patterning network.

Genotype References
wt [25], [28] hb 2 [25], [36] Kr 2 [25] pdm 2 [26], [41] cas 2 [26], [41] hb o.e. 1 [25] Kr o.e. [25], [40] pdm o.e. [26], [41] cas o.e. [26], [ Here i refers to one of each gene: i [ fhb, Kr, pdm, cas, xg. Here we have assumed that the expression noise comes from the transcription process (noise is incorporated only in the dynamics of {M i (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 {P i (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.
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 {P i } were discretized to 1 (0) for P i . P th (P i , P th ). The threshold P th was set as P th = 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.
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.
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 S pdm and S cas ), and the numbers of successful parameter sets are less than those obtained for the Drosophila network.
How is the robust nature of the Drosophila network implemented by these regulations? As seen above, the parameter values of S pdm and S cas (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 (J J pdm, Kr ,J J cas, pdm , and J J cas, hb ) and default promoter activities (S pdm and S cas ) (Fig. 8A-C). Figure 8A shows the fraction of successes for random assignments of parameter values under given strengths ofJ J pdm, Kr and S pdm . To score high reproducibility, S pdm must be large for smallJ J pdm, Kr , but need not to be large for sufficiently largẽ J J pdm, Kr . 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 betweenJ J cas, pdm and S cas (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 ofJ J cas, hb is small, S cas must be small to achieve a high fraction of successes (Fig. 8C). As D DJ J cas, hb D D becomes larger, a higher value of S cas 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  Table 4. doi:10.1371/journal.pcbi.1000760.g004 Table 3. Criterion for expression profile in each genotype.

Genotype
Criterion for the expression profile 1 wt 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].  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 (N d ). Here N d 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 N d = 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. doi:10.1371/journal.pcbi.1000760.g005  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. doi:10.1371/journal.pcbi.1000760.g006   Table 4. The temporal dynamics were tested for 50,000 trials. doi:10.1371/journal.pcbi.1000760.g008 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).
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 correla- tions 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.

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 t~t switch (0ƒt switch ƒt end). Once we assign the switching time of x expression (t switch ), 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 J hb, svp = 25, J hb, j = 0 (j = hb, Kr, pdm, cas, or x), and J k, svp = 0 (k = Kr, pdm, or cas) throughout this study. The time step at which we finish the simulation (t end ) was set as t end~1 2.
We thus investigated the behaviors of the remaining three factors (Kr, pdm, and cas) under the given regulatory interactions {J ij }. The total number of combinations of the parameters is 3 M |2 3 (the number of possible network architecture {J ij } multiplied by the number of default expression states fh i g), 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 fJ J ij g are continuous equivalents of {J ij } 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 {S i } are the default activities of the promoters. Transcription of a gene is induced when the total regulatory inputs become positive (S i z P jJ J ij P j w0), and is suppressed when they become negative (S i z P jJ J ij P j v0). In order to consider the effect of fluctuations on the expression dynamics, we introduced additive white Gaussian noise fj i (t)g: Sj i (t)j j (t 0 )T~s 2 i d ij d(t{t 0 ) (Eq (2)), where s i 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 (fJ J hb, i g~fJ J x, i g~0). To describe the expression change of hb and x, the promoter activities of these two are set as S hb .0 for tvt hb, off (S x .0 for tvt x, off ) and S hb = 0 for twt hb, off (S x = 0 for twt x, off ), respectively. The promoter activities of the others are always assumed to exist (S Kr , S pdm , and S cas .0). The noise intensities are also set as s i~s (.0) for tvt i, off and s i~0 for twt i, off (i = hb, x). Those of the other genes are s j~s (.0) (j = Kr, pdm, cas), Here we simply assume that the noise intensities of the genes take the same value s. The noise intensity s is set as s~0:05 in Figure 4, and s~0:08 in Figure 5. Noise intensity (horizontal axis) in Figure 6B means the value of s.

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 {S i } were assigned randomly within the defined ranges. The values of the matrixJ J ij È É were set to 0 when the corresponding regulations were absent (the corresponding element of the Boolean model takes J ij = 0) or assigned randomly when they are present (J ij =0). In order to confine our attention to the properties of network architectures, the other parameters (c M , c P , K M , and a) 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.
In the simulations, we found that the existence of positive selfregulation 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.