Inferring Drosophila gap gene regulatory network: a parameter sensitivity and perturbation analysis
 Yves FomekongNanfack^{1},
 Marten Postma^{1} and
 Jaap A Kaandorp^{1}Email author
DOI: 10.1186/17520509394
© FomekongNanfack et al; licensee BioMed Central Ltd. 2009
Received: 22 January 2009
Accepted: 21 September 2009
Published: 21 September 2009
Abstract
Background
Inverse modelling of gene regulatory networks (GRNs) capable of simulating continuous spatiotemporal biological processes requires accurate data and a good description of the system. If quantitative relations between genes cannot be extracted from direct measurements, an efficient method to estimate the unknown parameters is mandatory. A model that has been proposed to simulate spatiotemporal gene expression patterns is the connectionist model. This method describes the quantitative dynamics of a regulatory network in space. The model parameters are estimated by means of modelfitting algorithms. The gene interactions are identified without making any prior assumptions concerning the network connectivity. As a result, the inverse modelling might lead to multiple circuits showing the same quantitative behaviour and it is not possible to identify one optimal circuit. Consequently, it is important to address the quality of the circuits in terms of model robustness.
Results
Here we investigate the sensitivity and robustness of circuits obtained from reverse engineering a model capable of simulating measured gene expression patterns. As a case study we use the early gap gene segmentation mechanism in Drosophila melanogaster. We consider the limitations of the connectionist model used to describe GRN Inferred from spatiotemporal gene expression. We address the problem of circuit discrimination, where the selection criterion within the optimization technique is based of the least square minimization on the error between data and simulated results.
Conclusion
Parameter sensitivity analysis allows one to discriminate between circuits having significant parameter and qualitative differences but exhibiting the same quantitative pattern. Furthermore, we show that using a stochastic model derived from a deterministic solution, one can introduce fluctuations within the model to analyze the circuits' robustness. Ultimately, we show that there is a close relation between circuit sensitivity and robustness to fluctuation, and that circuit robustness is rather modular than global. The current study shows that reverse engineering of GRNs should not only focus on estimating parameters by minimizing the difference between observation and simulation but also on other model properties. Our study suggests that multiobjective optimization based on robustness and sensitivity analysis has to be considered.
Background
Gene regulatory networks (GRNs) play a fundamental role in body plan formation and development [1]. In most cases the published gene regulatory networks are based on experimental studies. Furthermore, various studies reveal that network dynamics depend on qualitative aspects (network structure) as well as the quantitative properties [2, 3]. An additional option is to study gene regulatory networks in mathematical models and simulation studies. Qualitative and quantitative models may provide insights in the causal relationships between components of a network as well as the mechanisms behind the network dynamics [2]. Qualitative models mainly focus on the inference of the wiring diagram and the understanding of the connection between the dynamics of the systems and the network structure [4, 5]. Understanding how the gene regulation leads to pattern formation requires the usage of quantitative models that allows inferring the network structure, but also simulating the spatiotemporal gene expression dynamics [6]. Many approaches have been used to model gene regulatory networks (see for a review [7]). In modelling processes from developmental biology especially the class of quantitative spatiotemporal, models of gene regulation is relevant. These type of models can potentially be linked with threedimensional biomechanical models of morphogenesis and provide new insights into developmental biology.
Although there is a large variety of formalisms used to describe quantitative model of GRN, the most commonly used is based on differential equations. This deterministic modelling framework allows describing spatiotemporal continuous system, where protein products are continuously produced through the balance of activation or repression. Models can be used in reverse engineering studies to infer the network structure and dynamics of the system. Quantitative spatiotemporal models of gene regulation are typically characterized by a large set of unknown parameters. Without a suitable method for reverse engineering the network and estimating the unknown parameter values these models have no practical use. Because of incomplete data, imprecise information concerning the interaction mechanisms and the large set of unknown parameter values; reverse engineering of the gene regulatory network is by far a straightforward task [8].
Reverse engineering of GRN capable of simulating spatiotemporal gene expression mainly consists in revealing the GRN structure that leads to the observed pattern. Optimization is therefore important [9, 10] and it is used to infer a gene network whether it is transcriptional or protein interaction. A system identification approach is used to select a model structure and by means of parameter estimation, the network topology is estimated from experimental data. The optimization problem consists in minimizing the difference between simulated expression profiles and available data in order to estimate the best circuit that predicts with very good accuracy the spatiotemporal gene expression. Based on the resulting parameterset, the network diagram is extracted and one tries to establish causalrelation of the dynamic mechanism that governs the pattern formation, using further analysis. The prediction is clearly linked to the quality and amount of the data; even with sufficient data, it is not guaranteed that the inferred network corresponds to the appropriate gene network that leads to the observed pattern [11, 12]. From a theoretical point of view, this question is a matter of the structural identifiability of the model [13]. Given a data set, is it possible to uniquely infer the network? Although there exist some theoretical methods [12] to investigate the identifiably before inferring the network, when confronted with a complex mechanism characterized by a multidimensional parameters space, the feasibility is still analytically complex. In many situations, one is confronted with an overfitting problem that leads to parameterset with very different quantitative values and even worse, yielding different network topologies having similar good realistic patterns.
To draw conclusions, one has to differentiate between circuits that are more likely to be biologically realistic.
It was shown that for certain experimental data that it is not possible to confirm whether the inferred model is really valid [14]. As long as the model is not contradicted by the given data, it is necessary to extend the validation test further. Especially when the model does not only predict a unique network, parameterset discrimination can be addressed using diverse models validation approaches. The quality of the circuits can for example be quantified by measuring the parameter reliability and sensitivity, the uniqueness of the predicted network, the model robustness and the predictability of the model. In this paper we focus on the model robustness against perturbation. In another paper (FomekongNanfack et al., in press.) we focus on the other aspects that can be used to quantify the quality of the circuits: methods for exploring the solution space, parameter correlation and asymptotical stability analysis.
By analyzing the robustness of the inferred network, we can test the quality of the circuits. The term robustness has various meanings, but in the current study, robustness is addressed as the ability of a system to maintain its mechanical and dynamical behaviour under perturbation. From a biological view, robustness is related to stability, homeostasis, canalization, redundancy, and plasticity [15–18], and can also be applied to dynamical process in development [19]. In simulations of dynamical developmental processes such as pattern formation, one should also expect robust behaviour within a certain extent in the model [20]. In the current paper, model robustness is addressed in two different perspectives. First, we investigate the quantitative robustness of the model towards internal fluctuations in expression level. It is known that presence of noise in gene regulation can lead to phenotype variation [21, 22]. There are some studies on the robustness of GRNs under the influence of molecular fluctuation [23, 24] and show the importance of noise and stochastic events [25]. In most cases, deterministic models are used to infer GRNs from noisy data. It has been shown in several studies that the robustness to noise mainly depends on the network structure instead of the parameter setting [26–28]. Therefore, one way to discriminate between circuits having different gene network but exhibiting the same pattern is to analyze their behaviour under noisy conditions. Modelsolutions showing more robustness or stability can be considered for further analysis. Second through simple parameter perturbation, we investigate how modelcircuits behave distinctly. This analysis allows us to identify the parameters that have the most significant influence on the model and to distinguish circuits that are less sensitive to overall perturbation. The combination of these two analyses allows one to discriminate between sets of circuits that are more robust, although they are obtained from the same model description and quantitative data.
ReverseEngineering the gap gene network
In the current paper, the case study early development along the anteroposterior (AP) of gap gene network of the Drosophila melanogaster is considered. In the gene regulatory network responsible for segmentation of the Drosophila melanogaster embryo [29–32] for most, if not all genes involved, experimental and bioinformatics studies are available. These studies not only give information about potential interactions between the different genes, but also provide spatiotemporal information about the gene expression patterns in the embryo. The segmentation genes leading to pattern formation, which occurs during the first stages of development after fertilization is controlled by a network of genetic interactions, with a cascade like modularity [33]. At the early stage, the maternal mRNA, located at the extremes of the egg will define the anteriorposterior axis of the embryo. Following fertilization, the process begins with the diffusion of maternal morphogen factors: Bicoid (Bcd), Nanos (Nos), Caudal (Cad) and maternal Hunchback (Hb), and the activation of the torso receptor on the poles. At cycle 13 (after 13 nuclear divisions), the gradient of these maternal gene products in the embryo promotes the transcription of the first zygotic genes during cycle 14: zygotic hb, giant (gt), Krüppel (Kr), knirps (kni) and tailless (tll); each forming expression domains in specific regions of the embryo [34]. This set of genes; the socalled gap genesare defined in broad domains along the anteriorposterior axis. Their transcription factors regulate the expression of another group of zygotic genes that comprises seven pairrule genes, which form a pattern of narrow stripes along the anteriorposterior axis. Finally, the pair rule genes regulate the formation of fourteen stripes by segment polarity genes. For review about the detailed mechanism, see [35–37].
For modelling the segmentation mechanism of the early Drosophila melanogaster embryo, two main formalisms have been proposed: logical formalism (tackling qualitative aspects) proposed by Sánchez and Thieffry [38], and continuous models proposed by Mjolsness et al. [39] used to obtain quantitative dynamics of a system. Following this formalism, Reinitz and coworkers formulated the problem as an inverse problem [40]. Given a complete mathematical model and sufficient accurate quantitative data, the parameters in the model can be estimated by optimization techniques, i.e., by fitting the model to the data. Except for box constraints, only little experimental information is used to constrain the parameter values in the model. This inverse modelling formulation of the problem leads to different gene circuits describing different aspect of the segmentation gene mechanism [41–44].
Using the gene circuit method proposed by Reinitz et al[40], Jaeger et al. [42, 45] have inferred a network by means of reverse engineering that can reproduce the measured spatial and temporal gene expression patterns. The gap gene model involves seven different genes, bcd, cad, hb, gt, Kr, kni and tll. The experimental data used to fit the model were obtained from the FlyEx database, where an extensive amount of accurate quantified spatiotemporal expression data for all genes is stored [46, 47]. A connectionist description was used to model the gene regulatory network. The number of parameters in this framework mounts up to 66 different unknowns for a network of six genes. To estimate these parameters the model was fitted to detailed spatiotemporal data using parallel simulated annealing (PLSA) [48–50], which is a technique that performs a global parameter search. Ten different gene circuits were obtained using this reverse engineering approach. Compared to the variance in the experimental data, all solutions fitted the dataset spatially and temporally accurately. Furthermore, the model reproduced the experimentally observed dynamic shift of the expression profiles and it was shown that diffusion is not the cause of the shift. Compared with experimental knowledge, for most interactions in the gene regulatory network the sign, i.e. activation, repression or nointeraction was reasonably well reproduced. Analysis of the 10 circuits suggested asymmetric repression between the gap genes, which was hypothesized to be the cause of the observed domain shift. This result has also been obtained by Perkins et al. [43] where parameter estimation was performed using a threestep strategy. Following [42], FomekongNanfack et al. [51] employed a different optimization technique to find parameter sets that accurately fitted the experimental dataset and lead to the similar networks as proposed by Jaeger et al. [45]. In search for a computational more efficient method, they developed a hybrid optimization algorithm [51]. Typically this approach first performs a global search based on evolutionary strategy followed by a local search. With this method, another 91 circuits were obtained that fitted the experimental dataset accurately. Recently, Ashyraliyev et al. [52] showed that the circuits obtained by FomekongNanfack et al. can be further improved using LevenbergMaquardt.
The gene circuit method previously used [40, 42, 43, 51] does not make any assumption about the network structure. All genetogene interactions are assumed to be plausible and minimizing the difference between observation and simulation drives the inference. The criterion for selecting acceptable circuits was based on a low root mean square error (RMS) and visual inspection of the simulated gene expression profiles (no major pattern defects noticeable). In the current case, a large set of circuits can be used to simulate the AP patterning of the Drosophila melanogater. Using a local sensitivity analysis, determinability of some parameters were studied in [52], where it was suggested that one could confirm on the nature of some biological interactions for parameters that were shown to be identifiable. Using a large set of parameters as base value, it is possible to examine the robustness of the model. Zak et al. [23] showed how input perturbations and stochastic gene expression influence the identifiability of a specific regulatory network given gene expression profiles and prior structural knowledge. Using as starting point the parametersset obtained in [42, 51] as base value, we discuss how does input perturbation and stochastic simulation allow a modelbased robustness analysis of a model that leads to multiple circuits.
Results
Gap gene circuits
Robustness towards fluctuations
Gene transcription, degradation and diffusion inherently are stochastic processes, which lead to fluctuations in gene expression levels. If expression levels are high, these fluctuations generally do not affect the time evolution of the system. In these cases the deterministic model and stochastic model will yield a similar result. However if the system is nonlinear and the fluctuations occur at an early time in development, where levels are still low, fluctuations may lead to different patterns. From a biological point of view this may not be favorable because gene pattern defects can lead to aberrations in development, hence development should be robust towards modest transient disturbances [53, 54]. Here we investigate if the circuits obtained from the optimization, which are based on a deterministic model, are robust towards fluctuations.
Analysis of fluctuations in the bicoid gradient [55] suggests that the expression level of bicoid is about five fold higher than the fluorescence unit. In our analysis we assume that the expression level of all genes is five fold higher than the fluorescence unit. We have used the Gillespie algorithm [56] to introduce fluctuations in the gene expression levels.
For each circuit we performed 100 stochastic simulations (see methods) and analyzed the quality of the final expression pattern at T = 68.1 min. In each run, by comparing the deterministic model with the stochastic model using a RMS cutoff criterion, we counted how many of domains cad, hbanterior (hba), hbposterior (hbp), Kr, gtanterior (gta), gtposterior (gtp), kni and tll were considered to be correct. For each circuit we also calculated an overall score by counting the runs where all expression domains formed correctly [see Tab.1 in Additional file 1, and Additional files 2, 3 w]. On average the domains Cad (99%) and Hbanterior (93%) show the best scores. The domains of kni (86%), Kr (80%) and gtposterior (81%) have an intermediate score. The domains gtanterior (66%), hbposterior (70%) and particularly tll (47%) have a low score. Most circuits have a very low overall score; only the top 25% of the circuits have an overall score higher than 35%. Most of the circuits have a score lower than 20%.
If all domains except one develop correctly the overall score is determined by that domain score, however if two or more independent domains have scores lower than 100%, the overall score will be lower than the individual domain scores. For example, if all domains would have a score of 99% and are independent the overall score would be 0.99^{8} * 100% = 92%. However we also observe interactions between domains, which leads to a higher overall score than would be expected if they behave independently. In these cases, if one domain yields a low score for a particular run, it is likely that another domain also yields a lower score for the same run. This is typically observed between adjacent domains. Notable examples are Kr with anterior gt, anterior hb with anterior gt, kni with posterior gt, and posterior gt with posterior hb and tll.
Additional file 4:stochastic simulation of circuit number 5 run 52. The movie shows the spatiotemporal simulation of circuit number 5 for the 52 stochastic run. (MP4 5 MB)
Additional file 5:stochastic simulation of circuit number 5 run 83. The movie shows the spatiotemporal simulation of circuit number 5 for the 83 stochastic run. (MP4 5 MB)
Additional file 6:stochastic simulation of circuit number 11 run 64. The movie shows the spatiotemporal simulation of circuit number 11 for the 64 stochastic run. (MP4 5 MB)
Additional file 7:stochastic simulation of circuit number 11 run 98. The movie shows the spatiotemporal simulation of circuit number 11 for the 98 stochastic run. (MP4 5 MB)
In Fig. 3A, B the final pattern formed during the stochastic simulations are shown for circuit 48 and 11. Circuit 48 has an overall pattern score of 62% and nr 11 has an overall score of 34%. In circuit nr 11 all domains except Kr (44%) and anterior gt (35%) score 100%. Kr and anterior gt show a strong interaction, if anterior gt disappears then Kr expands into the anterior gt domain. In circuit 48 this interaction between Kr and Gt at the anterior of the embryo is not present. In this circuit especially tll is not well defined. Fig. 4AB shows the trajectory of Kr and gt at nucleus 35, for circuit 48 in almost all runs gt and Kr develop correctly, however in circuit 11 there are two possible outcomes of the stochastic run. Here gt and Kr can develop correctly, however gt may also disappear and completely repressed by Kr. In this circuit there are two pathways for the system to evolve, which may lead to two stable points or a single stable point but two pathways. Jaeger et al. [45] suggested that nonoverlapping gap gene are mutually exclusive and have mutual repression. This result was confirmed for gt and Kr by Ashyraliyev et al. [52]. Circuit 48 and 11 both show this strong mutual repression but we believe that the bad stochastic score of anterior gt and Kr might be caused by the weak repression of gt by Tll.
In Fig. 3C, D the final pattern formed during the stochastic simulations are shown for circuit 2 and 41. Both circuits have a very low overall score of 12% and 0% respectively. In circuit 2, a very low score for hbp, gtp and tll mainly causes this and for circuit 41 all domains except cad are not well defined. In Fig. 4C, D the trajectories at nucleus 35 of hb and gt are shown. For circuit 2 the pathway is well defined, in most runs both gt and hb evolve similar to the deterministic model. However in circuit 41 the pathways are not well defined at all and show variability. In circuit 2, we see that Hb represses gt while Gt weakly activates hb, and inversely for circuit 41. It is suggested that overlapping gap genes do not activate each other [45, 52]. The poor robustness of these 2 circuits is therefore a consequence of the wrong connectivity of these 2 interactions.
In Fig. 3E, F the final pattern formed during the stochastic simulations are shown for circuit 5 and 24. Circuit 5 has a very low score of 10% and circuit 24 a score of 40%. Although circuit 5 has a low overall score both the kni and posterior gt domain score 100%, circuit 24 has a lower score for these domains. In Fig. 4E, F the trajectories are shown for kni and gt at nucleus 69. In circuit 5 the final points of the stochastic runs are very close to the final point of the deterministic run. Although the final points are well defined, the pathway, which is a shift of both domains, to these points is quite variable. In circuit 24 the same phenomenon is observed, only here Kni in some cases completely suppresses gt. Both circuits show gt activation by Kni (weak), which seems to be a wrong interaction [45, 52]. Circuit 5 shows repression of hb by Gt and weak activation of gt by Hb while circuit 24 shows the inverse mechanism. Based on the identifiability analysis obtained by Ashyraliyev et al. [52], it is suggested that Gt does not repress hb and Hb does not regulate gt. However, the determinability of these parameters has a very poor confidence and qualitative conclusion on these interactions is still ambiguous.
In Fig. 3G, H the final pattern formed during the stochastic simulations are shown for circuit 20 and 79. Circuit 20 has an overall score of 50% and circuit 79 has a very low score of 2%. In circuit 20 all domains except Kr and anterior gt are well defined. In circuit 79 all domains except cad are not well defined. In Fig. 4G, H the trajectory of hb and tll at nucleus 92 are shown. In circuit 20 most final points are very close to the deterministic model, however in circuit 79 almost all points are far away from the deterministic model. In this circuit tll domain is repressed by Hb and completely disappears and the hb domain continuous to grow. In some runs in circuit 20 (blue trajectories) we see the same tendency of continuous tll repression combined wit hb increase. Circuit 79 shows very inconsistent regulatory mechanism with respect to available literature [31, 32, 57, 58]. gt and hb show mutual activation, and Kni activates kr. These interactions seem to be the wrong regulatory mechanism, leading to stochastic instability and a very low pattern score.
Correlation between robustness and parameters
Robustness towards parameter perturbation
If a parameter in one of the circuits (considered as a model input) is perturbed by a certain amount the formation of the gene expression pattern (considered as the model output) will be perturbed as well. The amount by which a parameter can be increased and decreased before a certain pattern error (RMS) is reached can be used as a measure for pattern sensitivity with respect to that parameter. In order to get a better understanding of the model's sensitivity, for each circuit and each parameter in that circuit the sensitivity interval (SI) was calculated (see methods). From a biological point view pattern formation should be robust towards small perturbations in the parameters [59].
Circuit sensitivity
The standard deviation of the circuit parameters can be quite large compared to the average parameter value shown in Fig. 2. And furthermore, average parameter value and standard deviation strongly correlate (r = 0.81). In Fig. 9 the average log SIs versus standard deviation is plotted. Amongst similar parameter types (weights and nonweights) no correlation is observed between SI and standard deviation. This suggests that if a parameter is very variable across the circuits, i.e. it is not welldefined [60], it is not necessary less sensitive. A notable example is and , which both have a similar standard deviation but differ in sensitivity by 2 log units. Furthermore, we also observe that the standard deviation on average is much higher than the sensitivity interval.
Model sensitivity vs. pattern robustness
Discussion
The 101 gap gene circuits were previously obtained with inverse modelling [42, 51] of a seven gene (bcd, cad, hb, Kr, gt, kni and tll) connectionist model using detailed spatiotemporal gene expression data [61]. All circuits were selected only on the basis of a low RMS value (RMS ≤ 12), and in all cases the circuits were able to reproduce the dynamics and the patterns accurately. However, the parameter values of the circuits appeared not to be well defined, for some parameters multiple clusters were found, which represented multiple circuit topologies. Furthermore many parameters showed a single cluster with a high degree of scattering around the mean as shown in Fig. 2. Before any further tests, this result insinuates that for some parameters such as the nonregulatory parameters (showing major scattering), the notion of a global optimal parameter set makes no biological sense as suggest by Von Dassow et al. [27]. Also, the scatterings of most of the regulatory parameters are almost in a precise region of the search space, describing a precise qualitative behavior that can be linked to a potential interaction (repression, nointeraction or activation). The model always performs well and leads to very good pattern suggesting that the model is robust to the network topology as well as variation in individual parameters. Robustness to network topology is a behavior already observed by Von Dassow et al. [27] in the case of the segment polarity of Drosophila melanogaster.
To further analyze the properties of the model and to classify the quality of the circuits we conducted a perturbation analysis. We analyzed the robustness of 101 Drosophila melanogaster gap gene circuits using two different approaches, first the robustness of pattern formation towards intrinsic fluctuations in gene expression levels was investigated using a stochastic model and secondly the sensitivity of the circuits with respect to the parameters in the model was investigated using a simple parameter perturbation technique.
Robustness towards fluctuations
Robustness of the regulatory network is expected to be a fundamental property of the network, pattern formation should therefore not be compromised by small intrinsic fluctuations in gene expression levels. Therefore, the patterns obtained from stochastic simulations should statistically be similar to the pattern obtained from the deterministic simulation for robust circuits.
Robustness of domain formation
Looking at the evolution of the system in phase space, we found that the pathways in many cases were not well defined, which appeared to be linked to the existence of multiple attractors. The fluctuations allowed the system to jump from one attractor to another, causing the system to evolve into a pattern very different compared to the deterministic model. This suggests that by using the reverse engineering approach many circuits can be found, which are well defined for the deterministic case but not for the stochastic case. Multiple attractors can easily occur in a nonlinear model with many parameters, where overfitting may lead to many connections in the circuits. Also the type of feedback loops may be crucial, e.g. cad has a negative autoregulation, which reduces the effect of fluctuations, however all other genes have strong positive autoactivation [63, 64]; in these cases fluctuations are amplified if there are no negative feedback loops associated with the increase of the expression level.
The effect of certain parameters on robustness
The connectionist model contains a promoter threshold for each gene, which determines if gene production is on or off when there are no other inputs. In the current gap gene model the production of hb, Kr, gt and kni can only be turned on by maternal inputs or by autoactivation, therefore the promoter threshold are set to a fixed negative value, which was either H 3.5 or H = 2.5. Looking at the robustness of individual circuits, we find that circuits with promoter threshold set to H = 2.5 are more robust towards fluctuations compared to circuits with promoter threshold set to H = 3.5, and also their parameters are on average less sensitive. In general we find that weaker maternal inputs from bicoid and caudal increase the robustness with respect to fluctuations. Furthermore higher promoter rates, decay rates and diffusion coefficients improve the robustness towards fluctuations.
We have seen that circuits are relatively less sensitive with respect to diffusion coefficients, and with some extent, to decay time constant. Regarding the diffusion coefficient, this result confirms observations made elsewhere where it was suggested that the diffusion is not essential to explain precise gene expression pattern formation [42, 65], but the diffusion term does account for the boundary nuclei effect [66]. It was also shown by NussleinVolhard et al. [67] that the effect of diffusion is reduced with the exponential increase of the number of nuclei. Also, Gregor et al. [68] showed that diffusion coefficient does not scale with varying embryos length. The rest of the parameters have a mixed sensitivity behavior and do not fall into a specific category [see Tab. 2 in Additional file 1]. Some of the weights are very sensitive and some are not, idem for the production rate. Although we do not have precise biological explanations regarding the difference of the sensitivity behavior of the parameters, it might be interested to experimentally investigate if such difference is a property of the connectionist model or a characteristic of the regulatory mechanism. The sensitivity information can then guide the selection of the optimal mutation targets and thereby reduce the experimental effort. This validation could be done for example by measuring mRNA degradation rate of zygotic Hb in embryos with over expressed maternal hb, or by measuring the binding affinity in mutants. Also one could consider inducing genetic mutation to control kinetic parameters that can be measured [69].
We also find that certain parameters that regulate posterior tll and also some parameters that regulate anterior gt and Kr affect robustness of the corresponding domains. Furthermore we find that robustness and parameter sensitivity are linked. Especially for the group with a lower promoter threshold we find that circuits with a lower average sensitivity are more robust towards fluctuations. These results show that in the current gap gene model not just the network topology but to a large extent the precise value of the parameters determines robustness. We noticed that inputs of cad to other genes correlate negatively with robustness, these parameters are amongst the most sensitive in the model, with a  log SI in the order of 4. If we assume that fluctuations are proportional to the square root of the concentration level (this is a reasonable assumption for steady production) the relative fluctuation level is about , where n is the number of molecules. For n = 100, n = 1000 and n = 10000 the relative fluctuation level then is 10%, 3.2% and 1%. The most sensitive parameters in the model have a relative sensitivity of about 0.5  2% [see Tab. 2, Additional file 1]. This holds assuming that 20% of the circuit RMS is a critical cutoff.
Possible reasons for weak robustness
We conclude that the current connectionist model describing the gap gene segmentation obtained from reverse engineering techniques is not robust. Although some circuits appear more robust than others, most of the circuits are extremely sensitive towards gene expression level fluctuations and also pattern formation is very sensitive with respect to perturbation of a large number of parameters. In most cases the parameters causing this are linked to: Tll regulation, in some cases are linked to anterior Gt and Kr regulation and finally the promoter threshold (see Fig. 5). There may be a multiple reasons for weak robustness. First, the incompleteness of the model may be the cause. In the real system it is known that the terminal pathway with Tll, huckebein [70] and Torso also regulate gap genes, the latter two are missing in the current model. Furthermore, hb regulation is not only zygotic but also has a maternal component at the anterior end of the embryo that is also regulated by nanos, nanos and a maternal description hb mRNA are missing in the model. Furthermore, only a part of the embryo is considered, which may lead to boundary effects at the anterior end. Secondly, the current reverse engineering approach using a simple RMS optimization may lead to sensitive circuits, not only some time points are missing in the data also features in the data that are not represented in the model may be over fitted. Finally, the rather phenomenological approach of the connectionist model may also be a source of sensitivity because promoter threshold has a profound effect on robustness.
How to improve the network inference
In modelling processes from developmental biology especially the class of quantitative spatiotemporal models of gene regulation is relevant. This type of models can potentially linked with threedimensional biomechanical models of morphogenesis and provides new insights into developmental biology. Especially the quantitative spatiotemporal models of gene regulation are characterized by a large number of unknown parameters and an (infinite) class of potential solutions. So far, very few modelbased analysis method have been proposed to validate or invalidate models, especially for nonlinear and spatially distributed models [6, 71, 72].
The analysis presented in this paper shows the effect of missing genes, suggesting that it might be necessary to include more genes known to be involved in gap gene segmentation. Although data might not be available, one could consider combining synthetic data with observation. The hierarchical modular structure of the system should be considered within the model description [27]. Phenomenological models such as the connectionist model can generate realistic pattern, but they still fail at establishing the clear role of genetic perturbation. These models can be improved by considering using Hill type functions [27] and/or more detailed rate equations [65]
Another improvement would be to infer the network using stochastic models that translate the fluctuations in the system. It is computationally expensive to numerically solve such systems and it will require efficient optimizations methods [73].
To prevent overfitting, one might consider reducing the number of free parameters by adding more constraints known from previous experiments. In general the inverse problem that is only based on minimizing the RMS leads to overfitting, allowing for many solutions that all fit the dataset equally well but may not show correct behavior or properties beyond the dataset. It is therefore difficult to know, which solution is most similar to the real system. The different postoptimization analyses presented in this paper and in an accompanying manuscript ([74]) revealed some simple but efficient invalidation tests. Ideally, we would like to funnel all circuits into a pipeline of tests where the circuits are filtered and those passing all the tests are the robust and stable solutions. It may however turn out that none of the solutions will pass all the tests. Therefore we propose not only to minimize the RMS but also to incorporate other objectives into the optimization strategy. By introducing multiple objectives the solutions obtained should possess better properties and also have better defined parameters. The possible additional objectives for gene regulatory networks are:

Sensitivity constraints should lead to networks that are robust towards parameter or input perturbation. This could simply be addressed by incorporating a sensitivity analysis within the optimization procedure. Using the Levenbergmarquardt method for parameter estimation, Ashryalyev et al. [52] have investigated simultaneously while estimating the parameters, their identifiability. RodriguezFernandez et al. [75] have presented a hybrid method that can handle the sensitivity analysis while optimization. Using such method, we suggest attributing to each circuit a sensitivity value that determines the overall sensitivity of its parameter.

Noise robustness constraint is another improvement that should lead to circuits robust towards fluctuations (internal or external). This improvement can be achieved by inferring the network using stochastic models that include the fluctuations in the system. It is computationally expensive to numerically solve such systems and it will require efficient optimizations methods [73, 76]. Furthermore, It is known that patterning is insensitive to external fluctuation such as maternal gene expression. Consequently, an efficient model should also be able to simulate gapgene expression given noisy external Bcd expression. An efficient way to distinguish all the solutions obtained from the parameter estimation would then be to change the Bcd dosage during the optimization It was shown by Manu et al. [77] that Bcd variation is not the cause of the precise canalization of the gap genes but its a consequence of the cross regulation between zygotic gap genes [29]. In an optimization procedure, good circuits should be able to reproduce the patterns as well as the shift without showing any major fluctuation of the gap gene expression while varying the Bcd concentration.

Reducing network connectivity In our accompanying manuscript [74], we show that the long term dynamic revealed that in the presence of certain motifs, the circuits pattern converge to oscillatory attractors. This behavior is not desirable as we expect the patterns to converge to a stationary point being the steady state. Therefore, we would ideally ignore or penalized the circuits having this behavior. Also, we have shown that circuits with realistic topologies do not have these motifs and converge to the desired attractors. This could be obtained by either running the circuits beyond the data set to force asymptotic stability [66], or by adding an entropy function to prioritize circuits with the minimal connectivity. It was shown by Isalan et al. [78] that circuits where reduced connectivity are more robust to parameter variation in comparison to circuits with less connected members.
These different additional objectives would be integrated in a multiobjective optimization framework, which is an extension of the classic optimization of a singleobjective function (see Handl et al. for interesting review on Multiobjective optimization [79]). Many multiobjective algorithms exist [80, 81] such as method where multi objectives are transformed into one single objective [82–84], Pareto [85, 86] and nonPareto dominance approaches. To our knowledge, very few have used Multiobjective optimization for GRN inference [87, 88]. Based on recent reviews [9, 10], stochastic ranking evolutionary strategy (SRES) is one of the best methods for parameter estimation of biological problems. This method has been used to obtain some of the circuits analyzed in this article and obtained elsewhere [51]. The implementation of the SRES is suitable for multipenalties where the weighted sum of penalties is used [9] as a multiobjective constraints. The strength of the ES is its intrinsic parallel nature [89]. One alternative would be to use an island based ES where the number of island is determined by the different objectives. Each island minimizes the multiobjectives, but the distribution of the weight is different from island to island. Migration between islands will spread individuals with a particular strong property in the other islands. The main objective being the least square difference between the simulation and the data, it might be more practical to start with one population where the penalties have equal weights and later on, switch to islands once the LSE is acceptable. This would guarantee that in all islands, the individuals have at least a good RMS.
Conclusion
In this paper, we have provided a robust analysis of a model used to infer the gap gene regulatory network of Drosophila melanogaster. The model has been extensively used elsewhere [42, 43, 51] to simulate and provide some explanations concerning the regulatory mechanism that leads to precise pattern formation. Unfortunately, many assumptions were based on a limited number of circuits obtained using simulated annealing and previous researchers assumed that the model was correct with a correct topology. In this article, we have shown that the mathematical model leads to different circuits all capable of reproducing the quantitative spatiotemporal gene expression pattern. Consequently, it is difficult to decide based solely on the architecture which circuit is the correct one.
Robustness towards fluctuation has revealed that the overall gap gene domain tends to be poorly resistant to perturbation and this weak property could be related to some particular interactions predicted by the circuits. Furthermore, parameter perturbation analysis has shown that the circuits with lower sensitivity do not necessarily yield to robustness to fluctuation. The reason for these few exceptions is related to the promotor threshold and to local domain robustness, which can considerably affect the overall global robustness. Overall, the analysis shows that the network possesses modular robustness and some local properties may affect the robustness of a gene expression locally as show in Fig. 11. This feature is strongly related to the network topology as well as the interaction weights.
From a biological point of view, this paper has shown that it is difficult to relate the connectionist model with biological evidence. The model ability to simulate the gene expression does not necessarily provide meaningful information since alternative networks are predicted. Therefore, it would be interesting to test some of the different alternatives, especially when there is not yet any experimental evidence to invalidate. However, it was recently demonstrated that it is still possible to draw qualitative conclusions on the regulatory topology of the gap gene network [52]. The overall analysis has shown that based on the robustness toward gene expression fluctuations and parameter perturbations, it is possible to identify robust circuits as well as the parameters that are identifiable according to their sensitivity intervals. For a computational/system biologist, this shows that it is essential to further analyze a model prediction, when results are obtained from reverse engineering based on parameter estimation, since some of its properties may or may not invalidate the results. We have also provided some preliminary suggestions to efficiently improve the GRN inference to avoid reverse engineering that leads to circuits with different topologies, by controlling the optimization by means of multiobjectives minimization.
Methods
Inference of the Gap Gene model
The gap gene circuits analyzed in this paper were presented by Jaeger et al. [42] and FomekongNanfack et al. [51]. In both cases, the inference was performed using the same quantitative data, the same model description but different parameter estimation methods.
Quantitative data used are available online in the FlyEx database http://urchin.spbcas.ru/flyex/ or http://flyex.ams.sunysb.edu/flyex. The database presents a collection of spatiotemporal gene expression data obtained from fluorescently stained wildtype embryos for Eve protein and two other genes. Data were obtained by applying different image processing strategies [46]. The embryos are for different time ranging from cycle 7 to cycle 14A. In the simulation, data obtained at cycle 12 were used as initial conditions. For the genes Kr, gt, kni, Tll these are very close to zero and set to 0 in the simulations.
where N_{ g }denotes the number of genes or gene products involved and Φ is a sigmoid function with range (0,1). (t) represents the concentration level at time t of gene a in nucleus i with 1 ≤ i ≤ N and N the number of nuclei during a cleavage cycle. The concentration, , of the maternal gene bicoid is taken from experimental observations and is kept constant in time during the simulation. The parameters are: the regulatory weight matrix , describing the influence of gene b on gene a, the production rate R_{ a }, the activation threshold h_{ a }for Φ, the decay rate λ_{ a }, the diffusion coefficient D_{ a }, and the regulatory influence bcd_{ a }.
where E(θ) is given by Equation (2) and N_{ d }is the number of data points.
Sensitivity analysis and confidence interval
Stochastic behavior of the model
We implemented a nonlinear stochastic deterministic differential equation model of the gap gene where the dynamic is driven by Gillespie stochastic simulation algorithm [56]. The master equation is based on the differential equation describing the concentration change of a gene expression in a particular nucleus. Each equation describes three different reactions: protein production, protein decay and diffusion. Only one reaction can occur at a specific time. The Gillespie algorithm simulates the system by choosing first in a probabilistic manner which reaction occurs, and then estimates when does the next reaction will be realized. To speed the process, we have used a method [94] that only calculates and updates reactions that have changed. The probabilistic nature of the algorithm imposed us to run 100 simulations for each of the 101 gap gene circuits. All simulations used the same initial condition and the number of molecules is 5 times the deterministic case (based on [55]). For simplicity, analyses are only made on the final pattern (gastrulation time) by comparing the deterministic simulation with the stochastic one.
Declarations
Acknowledgements
This work was supported by the Netherlands Organization for Scientific Research, project NWOCLS 635.100.010 (Simulation of developmental regulatory networks", http://www.science.uva.nl/research/scs/3DRegNet) and by the EC (MORPHEX, NEST Contract No 043322). We used data from the FlyEx database http://flyex.ams.sunysb.edu. We thank Dr. Johannes Jaeger for his wise suggestions.
Authors’ Affiliations
References
 Davidson E: The regulatory genome  Gene regulatory networks in development and evolution. 2006, London: Academic PressGoogle Scholar
 Kitano H: Looking beyond the details: a rise in systemoriented approaches in genetics and molecular biology. Curr Genet. 2002, 41: 110. 10.1007/s002940020285zView ArticlePubMedGoogle Scholar
 Guet CC, Elowitz MB, Hsing W, Leibler S: Combinatorial Synthesis of Genetic Networks. Science. 2002, 296: 14661470. 10.1126/science.1067407View ArticlePubMedGoogle Scholar
 Julien Gagneur GC: From molecular networks to qualitative cell behavior. FEBS letters. 2005, 579 (8): 18671871. 10.1016/j.febslet.2005.02.007View ArticlePubMedGoogle Scholar
 Batt G, Ropers D, de Jong H, Geiselmann J, Mateescu R, Page M, Schneider D: Validation of qualitative models of genetic regulatory networks by model checking: analysis of the nutritional stress response in Escherichia coli. Bioinformatics. 2005, 21 (Suppl 1): i1928. 10.1093/bioinformatics/bti1048View ArticlePubMedGoogle Scholar
 Reeves GT, Muratov CB, Schupbach T, Shvartsman SY: Quantitative models of developmental pattern formation. Dev Cell. 2006, 11 (3): 289300. 10.1016/j.devcel.2006.08.006View ArticlePubMedGoogle Scholar
 de Jong H: Modeling and Simulation of Genetic Regulatory Systems: A Literature Review. J Comput Biol. 2002, 9: 67103. 10.1089/10665270252833208View ArticlePubMedGoogle Scholar
 D'haeseleer P: Reconstructing Gene Network from Large Scale. PhD thesis. 2000, University of New MexicoGoogle Scholar
 Moles CG, Mendes P, Banga JR: Parameter Estimation in Biochemical Pathways: A Comparison of Global Optimization Methods. Genome Res. 2003, 13 (11): 24672474. 10.1101/gr.1262503PubMed CentralView ArticlePubMedGoogle Scholar
 Ashyraliyev M, FomekongNanfack Y, Kaandorp JA, Blom J: Systems biology: parameter estimation for biochemical models. FEBS Journal. 2008, 276 (4): 886902. 10.1111/j.17424658.2008.06844.x.View ArticleGoogle Scholar
 Krishnan A, Giuliani A, Tomita M: Indeterminacy of Reverse Engineering of Gene Regulatory Networks: The Curse of Gene Elasticity. PLoS ONE. 2007, 2 (6): e562 10.1371/journal.pone.0000562PubMed CentralView ArticlePubMedGoogle Scholar
 Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP: Universally sloppy parameter sensitivities in systems biology models. PLoS Comput Biol. 2007, 3 (10): 18711878. 10.1371/journal.pcbi.0030189View ArticlePubMedGoogle Scholar
 Jaqaman K, Danuser G: Linking data to models: data regression. Nat Rev Mol Cell Biol. 2006, 7 (11): 813819. 10.1038/nrm2030View ArticlePubMedGoogle Scholar
 Poolla K, Khargonekar P, Tikku A, Krause J, Nagpal K: A TimeDomain Approach to Model Validation. IEEE Trans Automat Contr. 1994, 39: 951959. 10.1109/9.284871.View ArticleGoogle Scholar
 Carlson J, Doyle J: Complexity and robustness. Proc Natl Acad Sci. 2002, 99: 25382545. 10.1073/pnas.012582499PubMed CentralView ArticlePubMedGoogle Scholar
 Kitano H: Towards a theory of biological robustness. Mol Syst Biol. 2007, 3: 137 10.1038/msb4100179PubMed CentralView ArticlePubMedGoogle Scholar
 Jen E: Stable or robust? What's the difference?. Complexity. 2003, 8: 1218. 10.1002/cplx.10077.View ArticleGoogle Scholar
 Wagner A: Robustness and evolvability: a paradox resolved. Proc R Soc Lond B Biol Sci. 2008, 275: 91100. 10.1098/rspb.2007.1137.View ArticleGoogle Scholar
 Stelling J, Gilles E, Doyle F: Robustness properties of circadian clock architectures. Proc Natl Acad Sci USA. 2004, 101: 1321013215. 10.1073/pnas.0401463101PubMed CentralView ArticlePubMedGoogle Scholar
 Stelling J, Sauer U, Szallasi Z, Doyle FJ, Doyle J: Robustness of Cellular Functions. Cell. 2004, 118 (6): 675685. 10.1016/j.cell.2004.09.008View ArticlePubMedGoogle Scholar
 Thattai M, van Oudenaarden A: Intrinsic noise in gene regulatory networks. Proc Natl Acad Sci. 2001, 98: 86148619. 10.1073/pnas.151588598PubMed CentralView ArticlePubMedGoogle Scholar
 Blake W, Kaern M, Cantor C, Collins J: Noise in eukaryotic gene expression. Nature. 2003, 422: 633637. 10.1038/nature01546View ArticlePubMedGoogle Scholar
 Daniel E, Zak GEG, Schwaber JS, Francis J, Doyle I: Importance of Input Perturbations and Stochastic Gene Expression in the Reverse Engineering of Genetic Regulatory Networks: Insights From an Identifiability Analysis of an In Silico Network. Genome Research. 2003Google Scholar
 Jin Y, Lindsey M: Stability analysis of genetic regulatory network with additive noises. BMC Genomics. 2008, 9 (Suppl 1): S21 10.1186/147121649S1S21PubMed CentralView ArticlePubMedGoogle Scholar
 Rao CV, Wolf DM, Arkin AP: Control, exploitation and tolerance of intracellular noise. Nature. 2002, 420: 231237. 10.1038/nature01258View ArticlePubMedGoogle Scholar
 Barkai N, Leibler S: Robustness in simple biochemical networks. Nature. 1997, 387: 913917. 10.1038/43199View ArticlePubMedGoogle Scholar
 George von Dassow EMM, Eli Meir , Odell GM: The segment polarity network is a robust developmental module. Nature. 2000, 406: 188192. 10.1038/35018085View ArticleGoogle Scholar
 Kastner J, Solomon J, Fraser S: Modeling a Hox gene network in silico using a stochastic simulation algorithm. Developmental Biology. 2002, 246: 122131. 10.1006/dbio.2002.0664View ArticlePubMedGoogle Scholar
 Jackle H, Tautz D, Schuh R, Seifert E, Lehmann R: Crossregulatory interactions among the gap genes of Drosophila. Nature. 1986, 324 (6098): 668670. 10.1038/324668a0.View ArticleGoogle Scholar
 Reinitz J, Levine M: Control of the initiation of homeotic gene expression by the gap genes giant and tailless in Drosophila. Dev Biol. 1990, 140: 5772. 10.1016/00121606(90)90053LView ArticlePubMedGoogle Scholar
 Capovilla M, Eldon ED, Pirrotta V: The giant gene of Drosophila encodes a bZIP DNAbinding protein that regulates the expression of other segmentation gap genes. Development. 1992, 114: 99112.PubMedGoogle Scholar
 DE C, MSG C, X W, A P, D P: A selforganizing system of repressor gradients establishes segmental complexity in Drosophila. Nature. 2003, 426: 849853. 10.1038/nature02189View ArticleGoogle Scholar
 Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P: Molecular Biology of the Cell, Garland. 2002, FourthGoogle Scholar
 Rolando RiveraPomar HJ: From gradients to stripes in Drosophila embryogenesis: filling in the gaps. Trends in Genetics. 1996, 12 (11): 478483. 10.1016/01689525(96)100445View ArticleGoogle Scholar
 Akam M: The molecular basis for metameric pattern in the Drosophila embryo. Development. 1987, 101: 122.PubMedGoogle Scholar
 Ingham PW: The molecular genetics of embryonic pattern formation in Drosophila. Nature. 1988, 335 (6185): 2534. 10.1038/335025a0View ArticlePubMedGoogle Scholar
 St Johnston D, NussleinVolhard C: The origin of pattern and polarity in the Drosophila embryo. Cell. 1992, 68 (2): 201219. 10.1016/00928674(92)90466PView ArticlePubMedGoogle Scholar
 Sánchez L, Thieffry D: A logical analysis of the gap gene system. Theor Biol. 2001, 211: 114141. 10.1006/jtbi.2001.2335.View ArticleGoogle Scholar
 Mjolsness E, Sharp DH, Reinitz J: A connectionist model of development. J Theor Biol. 1991, 152 (4): 429453. 10.1016/S00225193(05)803911View ArticlePubMedGoogle Scholar
 Reinitz J, Sharp DH: Mechanism of eve stripe formation. Mech Dev. 1995, 49 (12): 133158. 10.1016/09254773(94)00310JView ArticlePubMedGoogle Scholar
 Sharp DH, Reinitz J: Prediction of mutant expression patterns using gene circuits. Biosystems. 1998, 47 (12): 7990. 10.1016/S03032647(98)000148View ArticlePubMedGoogle Scholar
 Jaeger J, Surkova S, Blagov M, Janssens H, Kosman D, Kozlov KN, Myasnikova E, VanarioAlonso CE, Samsonova M, Sharp DH, Reinitz J: Dynamic control of positional information in the early Drosophila embryo. Nature. 2004, 430 (6997): 368371. 10.1038/nature02678View ArticlePubMedGoogle Scholar
 Perkins TJ, Jaeger J, Reinitz J, Glass L: Reverse engineering the gap gene network of Drosophila melanogaster. PLoS Comput Biol. 2006, 2 (5): e51 10.1371/journal.pcbi.0020051PubMed CentralView ArticlePubMedGoogle Scholar
 Janssens H, Hou S, Jaeger J, Kim AR, Myasnikova E, Sharp D, Reinitz J: Quantitative and predictive model of transcriptional control of the Drosophila melanogaster even skipped gene. Nat Genet. 2006, 38 (10): 11591165. 10.1038/ng1886View ArticlePubMedGoogle Scholar
 Jaeger J, Blagov M, Kosman D, Kozlov KN, Myasnikova E, Surkova S, VanarioAlonso CE, Samsonova M, Sharp DH, Reinitz J: Dynamical analysis of regulatory interactions in the gap gene system of Drosophila melanogaster. Genetics. 2004, 167 (4): 17211737. 10.1534/genetics.104.027334PubMed CentralView ArticlePubMedGoogle Scholar
 Myasnikova E, Samsonova A, Kozlov K, Samsonova M, Reinitz J: Registration of the expression patterns of Drosophila segmentation genes by two independent methods. Bioinformatics. 2001, 17: 312. 10.1093/bioinformatics/17.1.3View ArticlePubMedGoogle Scholar
 Poustelnikova E, Pisarev A, Blagov M, Samsonova M, Reinitz J: A database for management of gene expression data in situ. Bioinformatics. 2004, 20 (14): 22122221. 10.1093/bioinformatics/bth222View ArticlePubMedGoogle Scholar
 Chu KW, Deng Y, Reinitz J: Parallel simulated annealing by mixing of states. J Comput Phys. 1999, 148 (2): 646662. 10.1006/jcph.1998.6134.View ArticleGoogle Scholar
 Lam J, Delosme JM: An Efficient Simulated Annealing Schedule: Derivation. Tech Rep 8816. 1988, Electrical Engineering Department, Yale, New Haven, CTGoogle Scholar
 Lam J, Delosme JM: An Efficient Simulated Annealing Schedule: Implementation and Evaluation. Tech Rep 8817. 1988, Electrical Engineering Department, New Haven, CTGoogle Scholar
 FomekongNanfack Y, Kaandorp JA, Blom J: Efficient parameter estimation for spatiotemporal models of pattern formation: case study of Drosophila melanogaster. Bioinformatics. 2007, 23 (24): 33563363. 10.1093/bioinformatics/btm433View ArticlePubMedGoogle Scholar
 Ashyraliyev M, Jaeger J, Blom JG: Parameter estimation and determinability analysis applied to Drosophila gap gene circuits. BMC Systems Biology. 2008, 2 (83):
 Goodwin BC, Kauffman S, Murray JD: Is Morphogenesis an Intrinsically Robust Process?. Journal of Theoretical Biology. 1993, 163: 135144. 10.1006/jtbi.1993.1112View ArticlePubMedGoogle Scholar
 Azevedo RB, Leroi AM: A power law for cells. Proc Natl Acad Sci USA. 2001, 98 (10): 56995704. 10.1073/pnas.091485998PubMed CentralView ArticlePubMedGoogle Scholar
 Wu YF, Myasnikova E, Reinitz J: Master equation simulation analysis of immunostained Bicoid morphogen gradient. BMC Syst Biol. 2007, 1: 52 10.1186/17520509152PubMed CentralView ArticlePubMedGoogle Scholar
 Gillespie DT: Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry. 1977, 81 (22): 23402361. 10.1021/j100540a008.View ArticleGoogle Scholar
 Tautz D: Regulation of the Drosophila segmentation gene hunchback by two maternal morphogenetic centres. Nature. 1988, 332 (6161): 281284. 10.1038/332281a0View ArticlePubMedGoogle Scholar
 Eldon E, Pirrotta V: Interactions of the Drosophila gap gene giant with maternal and zygotic patternforming genes. Development. 1991, 111 (2): 367378.PubMedGoogle Scholar
 Umulis D, O'Connor MB, Othmer HG: Robustness of embryonic spatial patterning in Drosophila melanogaster. Curr Top Dev Biol. 2008, 81: 65111. full_textView ArticlePubMedGoogle Scholar
 Sun N, Sun NZ, Elimelech M, Ryan JN: Sensitivity analysis and parameter identifiability for colloid transport in geochemically heterogeneous porous media. Water Resources Research. 2001, 37 (2): 209222. 10.1029/2000WR900291.View ArticleGoogle Scholar
 Myasnikova EM, Samsonova AA, Samsonova MG, Reinitz J: Spatial registration of in situ gene expression data. Molecular Biology. 2001, 35 (6): 955960. 10.1023/A:1013215108374.View ArticleGoogle Scholar
 Bronner G, Jackle H: Control and function of terminal gap gene activity in the posterior pole region of the Drosophila embryo. Mech Dev. 1991, 35 (3): 205211. 10.1016/09254773(91)900193View ArticlePubMedGoogle Scholar
 Dunlap JC: Molecular Bases for Circadian Clocks. Cell. 1999, 96 (2): 271290. 10.1016/S00928674(00)805668View ArticlePubMedGoogle Scholar
 Cross FR, Siggia ED: Shake It, Don't Break It: Positive Feedback and the Evolution of Oscillator Design. Developmental Cell. 2005, 9 (3): 309310. 10.1016/j.devcel.2005.08.006View ArticlePubMedGoogle Scholar
 Alves F, Dilao R: Modeling segmental patterning in Drosophila: Maternal and gap genes. J Theor Biol. 2006, 241 (2): 342359. 10.1016/j.jtbi.2005.11.034View ArticlePubMedGoogle Scholar
 Gursky V, Kozlov K, Samsonov A, Reinitz J: Model with asymptotically stable dynamics for Drosophila gap gene network. Biophysics. 2008, 53 (2): 164176. 10.1134/S0006350908020085.View ArticleGoogle Scholar
 NussleinVolhard C: Gradients that organize embryo development. Sci Am. 1996, 275 (2): 5455.View ArticlePubMedGoogle Scholar
 Gregor T, Bialek W, de Ruyter van Steveninck RR, Tank DW, Wieschaus EF: Diffusion and scaling during early embryonic pattern formation. Proc Natl Acad Sci USA. 2005, 102 (51): 1840318407. 10.1073/pnas.0509483102PubMed CentralView ArticlePubMedGoogle Scholar
 Feng XJ, Hooshangi S, Chen D, Li G, Weiss R, Rabitz H: Optimizing genetic circuits by global sensitivity analysis. Biophys J. 2004, 87 (4): 21952202. 10.1529/biophysj.104.044131PubMed CentralView ArticlePubMedGoogle Scholar
 Weigel D, Jurgens G, Klingler M, Jackle H: Two gap genes mediate maternal terminal pattern information in Drosophila. Science. 1990, 248 (4954): 495498. 10.1126/science.2158673View ArticlePubMedGoogle Scholar
 ElSamad H, Prajna SAP, Doyle J, Khammash M: Advanced methods and algorithms for biological networks analysis. Proc IEEE. 2006, 94: 832853. 10.1109/JPROC.2006.871776.View ArticleGoogle Scholar
 Prajna S: Barrier certificates for nonlinear model validation. Automatica. 2006, 42: 117126. 10.1016/j.automatica.2005.08.007.View ArticleGoogle Scholar
 Tomshine J, Kaznessis YN: Optimization of a stochastically simulated gene network model via simulated annealing. Biophys J. 2006, 91 (9): 31963205. 10.1529/biophysj.106.083485PubMed CentralView ArticlePubMedGoogle Scholar
 FomekongNanfack Y, Postma M, Kaandorp JA: Inferring Drosophila gap gene regulatory network: pattern analysis of simulated gene expression profiles and stability analysis. [To appear in BMC System Biology].
 RodriguezFernandez M, Mendes P, Banga JR: A hybrid approach for efficient and robust parameter estimation in biochemical pathways. Biosystems. 2006, 83 (23): 248265. 10.1016/j.biosystems.2005.06.016View ArticlePubMedGoogle Scholar
 Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev Genet. 2009, 10 (2): 122133. 10.1038/nrg2509View ArticlePubMedGoogle Scholar
 Manu , Surkova S, Spirov AV, Gursky VV, Janssens H, Kim AR, Radulescu O, VanarioAlonso CE, Sharp DH, Samsonova M, Reinitz J: Canalization of Gene Expression in the Drosophila Blastoderm by Gap Gene Cross Regulation. PLoS Biology. 2009, 7 (3):
 Isalan M, Lemerle C, Serrano L: Engineering gene networks to emulate Drosophila embryonic pattern formation. PLoS Biol. 2005, 3 (3): e64 10.1371/journal.pbio.0030064PubMed CentralView ArticlePubMedGoogle Scholar
 Handl J, Kell DB, Knowles J: Multiobjective optimization in bioinformatics and computational biology. IEEE/ACM Trans Comput Biol Bioinform. 2007, 4 (2): 279292. 10.1109/TCBB.2007.070203View ArticlePubMedGoogle Scholar
 Deb K: MultiObjective Optimization using EvolutionaryAlgorithms. 2001, WileyGoogle Scholar
 Colette Y, Siarry P: Multiobjective Optimization, Principles and Case Studies. 2003, SpringerGoogle Scholar
 Syswerda G, Palmucci J: The application of genetic algorithms to resource scheduling. Fourth Int Conf on Genetic Algorithms ICGA'4. Edited by: Belew R, Booker L. 1991, 502508. San Mateo, California.: Morgan Kaufmann PubGoogle Scholar
 Jakob W, GorgesSchleuter M, Blume C: Application of genetic algorithms to task planning and learning. Parallel Problem Solving from Nature PPSN'2, LNCS. Edited by: Manner R, Manderick B. 1992, 291300. Amsterdam. NorthHollandGoogle Scholar
 Jones G, Brown R, Clark D, Willett P, Glen R: Searching databases of twodimensional and threedimensional chemical structures using genetic algorithms. Fifth Int. Conf. on Genetic Algorithms. Edited by: Forrest S. 1993, 597602. San Mateo, California.: Morgan KaufmannGoogle Scholar
 Goldberg DE: Genetic Algorithms in Search, Optimization, and Machine Learning. 1989, AddisonWesley ProfessionalGoogle Scholar
 Bentley P, Wakefield J: Soft Computing in Engineering Design and Manufacturing. 1997, 231240. chap. Finding acceptable Paretooptimal solutions using multiobjective genetic algorithmsGoogle Scholar
 Praveen K, Sanjoy D, Stephen W, L RJ: A multiobjective GAsimplex hybrid approach for gene regulatory network models. IEEE Proceedings of the 2004 congress on evolutionary computation: CEC 2004. 2004, 20842090.Google Scholar
 Liu PK, Wang FS: Inference of biochemical network models in Ssystem using multiobjective optimization approach. Bioinformatics. 2008, 24 (8): 10851092. 10.1093/bioinformatics/btn075View ArticlePubMedGoogle Scholar
 Lohmann R: Application of Evolution Strategy in Parallel Populations. Parallel Problem Solving from Nature  Proceedings of 1st Workshop, (PPSN 1), Volume 496 of Lecture Notes in Computer Science. Edited by: Schwefel HP, Mäanner R. 1991, 198208. Berlin, Germany: SpringerVerlagGoogle Scholar
 Foe VE, Alberts BM: Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in Drosophila embryogenesis. J Cell Sci. 1983, 61: 3170.PubMedGoogle Scholar
 Runarsson TP, Yao X: Search biases in constrained evolutionary optimization. IEEE Transactions on Systems, Man, and Cybernetics, Part C. 2005, 35 (2): 233243. 10.1109/TSMCC.2004.841906.View ArticleGoogle Scholar
 Lewis RM, Shepherd A, Torczon V: Implementing generating set search methods for linearly constrained minimization. Tech Rep WMCS200501. 2005, [Revised July 2006], Department of Computer Science, College of William & MaryGoogle Scholar
 Holland JH: Genetic Algorithms. Sci Am. 1992, 267: 6672.View ArticleGoogle Scholar
 Gibson M, Bruck J: Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels. Journal of Physical Chemistry A. 2000, 104 (9): 18761889. 10.1021/jp993732q.View ArticleGoogle Scholar
 Kosman D, Small S, Reinitz J: Rapid preparation of a panel of polyclonal antibodies to Drosophila segmentation proteins. Dev Genes Evol. 1998, 208: 290294. 10.1007/s004270050184View ArticlePubMedGoogle Scholar
 Myasnikova E, Kosman D, Reinitz J, Samsonova M: SpatioTemporal Registration of the Expression Patterns of Drosophila Segmentation Genes. Proceedings of the Seventh International Conference on Intelligent Systems for Molecular Biology. 1999, 195201. AAAI PressGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.