BACKGROUND
gene regulatory networks  are essential for controlling cellular metabolism and the organismal development.
under the command of transcription factors , each gene influences the activity of the cell by generating messenger rna  that guides the synthesis of proteins by ribosomes in the cytoplasm .
gene network modeling uses gene expression data to characterize the phenotypic behavior of a system under study.
with the reconstructed networks, biologists can generate and test hypotheses to further understand the complex phenomena that occur in nature systems and to explore the dynamics of those systems.
modeling grns manually on the basis of the experimentally measured time-series data takes a considerable amount of time.
therefore, an automated reverse-engineering procedure is recommended  <cit> .
this procedure involves altering the gene network in some way, observing the outcome, and using computational methods to infer the underlying principles of the network.
to derive a realistic model, available domain knowledge  can be integrated into the computational methods.
figure  <dig> illustrates the general procedure of inferring grns from quantitative expression data.
the inferring/modeling block of this figure indicates the computational procedure used to derive network parameters for a given model, to build and simulate the model, and to evaluate it by comparing the behavior of the inferred model with the original data set.
in addition to making use of expression data, one recently developed strategy combines information from various sources to narrow down the search space in the network.
this strategy shortens the time and effort required for the validation and discovery of networks.
for example, if some gene names are known then they can be mapped into knowledge bases  to extract biological knowledge  that can be used to determine the network structures.
our goal in this study is to establish a methodology for network inference and to investigate aspects of network dynamics that have not been addressed previously.
to infer a network with desired system behavior, the most important steps are to select a network model and then to fit the network's structural parameters to the available expression data.
many models have been proposed to address different levels of biological details, ranging from the very abstract  to the very concrete .
abstract models are easy to simulate and therefore less computationally taxing, but it has been proven that they are not able to capture certain system behaviors.
in contrast, concrete models are more suitable for simulating biochemical process realistically, but, because of their computational complexities, these models can only be applied to small systems.
to capture the underlying physical phenomena of a gene network we took one of the most popular and well-studied concrete models, the s-system model, as an example of a biological network, and we used both simulated and collected real gene expression data to reconstruct the model.
the s-system model is a type of ordinary differential equation  model.
it consists of a particular set of tightly coupled odes in which the component processes are characterized by power law functions  <cit> .
the system structure of an s-system is described by the following equation:
dxidt=αi ∏j=1nxjgi,j-βi ∏j=1nxjhi,j
here, xi is the expression level of gene i, and n is the number of genes in a genetic network.
the non-negative parameters αi and βi are rate constants that indicate the direction of mass flow.
the real-number exponents gi,j and hi,j are kinetic orders that reflect the strength of interactions from genes j to i. the above set of parameters defines an s-system model.
to infer an s-system model is, thus, to estimate all of the 2n  parameters simultaneously.
although non-linear odes can more accurately model the system dynamics of gene networks, they are difficult to solve by traditional optimization techniques  <cit> .
to reduce the complexity of such models, di bernardo et.
al developed an approach that involves taking a series of steady-state gene expression measurements following transcriptional perturbations and reducing the linear odes to a linear regression problem that can be solved with relevant techniques  <cit> .
in addition, some researchers employed the strategies of gene clustering for dimension reduction and gene classification for identification of expression patterns to reduce task complexity  <cit> .
as evolutionary algorithms  have been widely used to solve many difficult optimization tasks with good results, they have been suggested as a way to infer gene networks  <cit> .
among many studies on this topic, the most relevant are those involving the use of eas to infer s-system models .
in these studies, the network parameters of the s-system were arranged as a string of floating numbers  that could be evolved by genetic operators.
in recent studies of ea-based parameter inference, one critical problem, network robustness  <cit> , has not been addressed.
it is important to investigate the effect of network parameter perturbations on the overall system.
these parameters not only are the numerical components of the model, but also represent the activities or interactions among proteins, transcription factors, and mrnas in the gene network.
each network parameter has an important role in determining the behavior of a biological system.
parameters that are very sensitive to variation can introduce fragility into the system.
in fact, recent research has found that parameters can define the dynamics of a model, and more important, studies have shown that enforcing constraints on the parameters can limit  model dynamics  <cit> .
specifically, by measuring and analyzing the variations of network parameters and their effects on relevant genes, the gene-gene interactions of an inferred model can be interpreted.
to ensure the robustness of the inferred network and to further investigate gene interactions, it is important to derive an acceptable value range for each parameter and to restrict the parameter's value to the specified range during the network reconstruction process  <cit> .
therefore, we take parameter sensitivity into consideration in the network reconstruction procedure so that robust results can be obtained.
sensitivity analysis  is an indispensable technique that can help researchers to investigate the parameter robustness properties of an inferred network.
by varying the parameter values within a certain range and performing statistical calculations to measure the system instability , researchers can identify the critical parameters or discover seemingly unimportant parameters that may have a positive or negative influence on a network.
sa techniques have been used in many biological studies, such as genetic circuit design  <cit> , mammalian circadian clock modeling  <cit> , and target prediction in signaling pathways  <cit> .
in general, the sensitivity of a parameter is defined as  <cit> :
spm=∂m/m∂p/p=percentage change in mpercentage change in p
here, p represents the parameters that are varied in a given range, m is the mathematical function describing the system behavior, and ∂m means the change in m due to the value changed in ∂p with respect to p.
in a dynamic system, there are two types of methods often used to conduct sa: local sa and global sa.
local sa can measure the sensitivity of a parameter from a given range even if the system structure is unknown, but it only considers one parameter at a time and ignores the interactions between parameters  <cit> .
alternatively, global sa can examine parameter interactions with different parameter magnitudes simultaneously  <cit> , but this can only be done if the system structure is known in advance.
the above two methods have been compared extensively  <cit> .
it should be noted that neither local sa nor global sa is superior: the choice of method depends on the specific application and the information available.
in the biological cases that involve the use of sa techniques, the sensitivity of parameters should be investigated from a biological point of view in order to validate the results.
from the perspective of parameter dynamics, we propose a new approach to infer robust networks.
our approach includes two major parts.
the first part is a sensitivity analysis procedure that selects sensitive network parameters, determines value ranges for them, and then sends the parameters with their constraints to the inference mechanism.
the second part is an evolutionary method responsible for inferring networks from the results obtained in the first part.
this approach, to our knowledge, is the first work in network modeling that integrates sensitivity analysis into an inference algorithm to consider both internal  and external  characteristics of the inferred network.
to validate the proposed approach, a series of experiments has been conducted, in which the proposed sa procedure was coupled with different eas.
the results have been analyzed, and they show that this approach can be used to infer robust networks successfully from artificial and real gene expression profiles.
methods
parameter sensitivity analysis
in the process of inferring a gene network from expression data the genes interact with each other, and the network structure is generally unclear during the modeling process.
most studies of global sa use known pathways as examples and can therefore easily select the critical parameters and control the total number of parameters  for further analysis  <cit> .
in the type of network inference problem addressed here, it is not possible to identify and select the most important parameters in order to perform global sa, because no prior information about network structure is available.
to consider multiple network parameters simultaneously in such situations, we present a new method that is a modification of a widely used global sa technique, called multi-parameter sensitivity analysis   <cit> .
mpsa is based on monte carlo simulations and on quantitative comparisons of cumulative frequency via corresponding pearson correlation coefficients among parameters to identify those that are most sensitive.
it can be used to explore critical genetic reactions and to examine the robustness of a gene network as other global sa methods do.
in the method described here, the objective function of the sa method is the fitness function that is to be optimized in the evolutionary algorithm, which is the mean squared error over the time course:
∑i=1n∑t=1txia-xidxid <dig>
in this expression, xid is the desired expression level of gene i at time t, xia is the actual value obtained from the inferred model, n is the number of genes in the network, and t is the number of time points used to measure gene expression data during the period.
our modified method, m-mpsa, inherits qualities from mpsa but has several advantages.
first, mpsa is easy to implement and modify.
it allows users to select one parameter at a time or to consider multiple parameters at the same time to perform sensitivity analysis.
this characteristic is crucial because, when the network structure is unknown or uncertain, selecting one parameter at a time is a good way to reduce the complexity of the search space.
researchers can then concentrate on each parameter and temporarily neglect the combination sets at the selection phase.
second, compared with most of the global sa methods, the mathematical equations of mpsa are more manifest and succinct, and the parameter sensitivity is thus easy to calculate.
finally, the value ranges of parameters used in mpsa are the same as in an ode model.
this means that the sensitivity of each parameter can be directly mapped into an ode model without any extra computation.
in previous studies, researchers have shown that the most influential  parameters play key roles in modeling a grn; that is, a model's behavior varies largely because of its sensitive parameters  <cit> .
therefore, our current study draws on the characteristics of sensitive parameters, and it incorporates the sensitivity analysis method into the incremental evolution approach to infer a grn.
the proposed m-mpsa method is described in the sensitivity analysis algorithm below.
it uses an iterative process to calculate the sensitivity of each parameter, and it then ranks the sensitivities of all parameters.
the input of this algorithm is the set of network parameters to be determined; and the output, a list of parameters ranked by their sensitivity values.
with the parameter sensitivities ranked, the evolutionary algorithm can then infer robust solutions by exploiting the sensitive parameters that are likely to significantly influence the genetic model as a whole.
sensitivity analysis algorithm()
step 1: select one parameter, i, at a time.
.
step 2: set the parameter range ri for each s-system parameter, i. initially, the commonly used search regions  and  are taken as the initial settings for the parameters of kinetic orders and rate constants respectively, as suggested in  <cit> .
these settings are used to define the parameter range ri, by adding/subtracting an amount of one third of the initial settings to/from the current value of parameter i.
step 3: for each ri, a set of random values are created with a uniform distribution within the specified range.
each random value for parameter i, together with the other  parameters, constitutes a candidate solution  and calculate the objective function value for each of the random values generated in step  <dig>
step 5: determine whether the objective function value of each random point obtained in step  <dig> is acceptable or unacceptable by a given threshold cr .
the parameter value of a random point is classified as unacceptable if its corresponding objective function value is greater than cr; otherwise, the parameter value is acceptable.
go back to step  <dig> until all parameters have been dealt with.
step 6: calculate the sensitivity for parameter i by using its cumulative frequency .
cfi measures the similarity of two statistical distributions that are formed by the acceptable and unacceptable values produced from steps 3- <dig> for parameter i. the parameters are ranked according to their cf, and the ones with relatively low cf are considered to be sensitive.
step 7: output a parameter list in which the parameters are sorted by their cf values.
in addition to presenting the m-mpsa algorithm, we provide here a walkthrough example to illustrate how it operates in practice.
in the example, we randomly choose two parameters from one of the datasets  used in the experimental section, and we calculate the sensitivity values of these parameters.
the dataset consisted of ten nodes, and there were altogether  <dig> parameters to be estimated.
without losing generality, we assume that the algorithm selected the parameter p <dig> , and then set the range r <dig> for p <dig> .
as mentioned above, we used one third of the original search region  to set up the bounds for each parameter.
hence, if p <dig> had a value of  <dig> , then the lower bound for this parameter was - <dig>  )), and the upper bound was  <dig> .
then,  <dig> random values were generated within the range  for p <dig> , and their corresponding objective  function values were calculated .
next, a threshold cr was determined to categorize the above  <dig> evaluation results as acceptable or unacceptable.
if the current best result  was  <dig> , then a threshold of  <dig>   was used to determine whether a result was acceptable or not .
after that, two statistical distributions were constructed for the acceptable and unacceptable results identified above.
for statistical purposes, the value range  was divided into ten intervals, and for each interval the numbers of acceptable and unacceptable values were counted.
figure  <dig>  shows the distributions for p <dig>  in which acceptable and unacceptable numbers were marked.
the cumulative frequency distributions were then constructed accordingly, as presented in figure  <dig> .
finally, the sensitivity value of p <dig> was obtained by measuring the correlation coefficients of the two  cumulative frequency distributions .
in this case, the sensitivity value of p <dig> was  <dig> .
a similar procedure was performed for the other parameter, p <dig>  and a sensitivity value of  <dig>  was obtained.
because p <dig> had a sensitivity value lower than that of p <dig>  it was considered to be more sensitive than p <dig>  in this way, after sensitivity values of all parameters were calculated, they were ranked and sent to the incremental evolution algorithm.
as mentioned above, more sensitive parameters are more likely to have a significant influence on the entire genetic model.
therefore, if no prior knowledge can be applied to the reconstruction of a large grn, researchers can conduct sensitivity analysis for parameter selection to guide the network inference method, and observe how numerical quantities within different intervals can formulate candidate results .
taking p <dig> as an example, the acceptable values for this parameter fell in the range  that consisted of  <dig> acceptable network models  in total.
we can exploit this result to infer p <dig> and explore the relationships of this parameter with other genetic reactions .
in fact, researchers have shown that the most influential  parameters play important roles in modeling a grn.
that is, the behavior variation of a network model largely depends on its sensitive parameters  <cit> .
in figure  <dig>  we can see that the variation of the parameter value range influenced the system dynamics.
that is, the system robustness was significantly affected by the value range specified for a parameter.
for instance, the two distributions for p <dig> in figure  <dig>  show that, if the value range of p <dig> changed from  to , then the acceptable values for p <dig> changed from  <dig> to  <dig>  this means that the inference algorithm can maintain the system robustness if it can locate p <dig> within the interval .
the above analysis concludes that the proposed sa method can be used to find sensitive parameters and to derive suitable value ranges as search constraints, which can then be used in the inference algorithm to construct robust networks.
evolutionary algorithm for parameter optimization
to derive the network parameters of s-system model, we also implement an evolutionary algorithm  for parameter optimization to work with the above sa method.
ea is a population-based approach that evaluates many solutions simultaneously in the search space, and it is likely to find a global solution for a given problem.
recently, a new population-based optimization technique, particle swarm optimization , was proposed as an alternative to the traditional eas.
pso tries to mimic the goal-seeking behavior of biological swarms.
the standard pso algorithm contains a set of particles and operates in an iterative manner.
each particle is characterized by its position and velocity, and it moves in the solution space.
the position of each particle represents a potential solution that is evaluated by a predefined fitness function.
pso has some attractive characteristics.
in particular, it has memory, so that knowledge of a good solution can be retained by all particles .
during the iterative search process, each particle remembers its previous best position and the best position of any particle in the swarm.
then the particle uses the position information to modify its position and velocity, and continues its movement in the search space.
the details of the pso algorithm can be found in  <cit> .
some performance comparisons between pso and the most popular ea method, the genetic algorithm , have been made, underscoring the reliability and the convergence speed of both methods.
however, the result tends to be inconclusive: each of the algorithms has shown better performance than the other for some particular applications.
consequently, hybrid techniques were proposed to effectively exploit the qualities and the uniqueness of the two methods.
it is now commonly agreed that the hybrid methods can lead to further performance improvements  <cit> .
therefore, in this work, we develop a hybrid method to exploit the solution memory of pso and the genetic operations of ga for network inference.
as has been mentioned, our primary goal is to investigate how parameter sensitivity can be used to guide an ea and derive robust networks.
to concentrate on this issue, we have chosen to implement a popular ga-pso method  for parameter optimization, and we have built our sa method upon it.
breeding swarms  <cit>  is a frequently used ga-pso hybrid method that has demonstrated its good performance in many benchmark test functions.
therefore, we adopt this method for our network inference application.
the hybrid ga-pso method used in this work is illustrated in figure  <dig>  initially, a population is randomly generated and evaluated.
the individual solutions are ranked according to their fitness values, and then they are divided into two parts for running the pso and ga methods separately.
as shown in figure  <dig>  the best individuals % of the whole population) are preserved and enhanced by the pso procedure.
they are then sent to the next generation.
meanwhile, on the right-hand side of the figure, the remaining individuals  are discarded.
to replace the removed individuals, a tournament selection scheme is used to choose the same number of individuals from the best ones , and the selected individuals are used to create new individuals by the ga procedure.
in our implementation, we take a direct encoding scheme to represent solutions for both the ga and pso parts, in which the network parameters  are arranged as a linear string chromosome of floating-point numbers.
the goal here is to minimize the accumulated discrepancy between the gene expression data recorded in the dataset and the values produced by the inferred model.
therefore, the error  function defined previously is used directly for performance measurement.
for the pso part, the equations for updating the particle's velocity and position are the same as the ones listed in the original pso work  <cit> , and for the ga part, the operations of crossover and mutation described in  <cit>  are used.
sensitivity-based incremental evolution method for network inference
as described in the first section, to infer the s-system model is to determine the 2n  parameters simultaneously.
solving this high-dimensional problem is difficult, especially when the complexity of regulation increases along with the number of genes involved.
one promising approach to this problem is to adopt the concept of incremental evolution to infer the large number of parameters.
the underlying principle of incremental evolution  <cit>  is that a population is first evolved to solve an easier version t' of the original complex task t, in which the solution region of t is more accessible from region t'.
more task versions with incremental complexity can be arranged so that the original task can be achieved progressively.
evidently, the main task involved in implementing incremental evolution is the formulation of a scheme to transfer the goal task into another more evolvable task.
in the process of task transformation, the structure of the environment and the overall goal must be preserved.
this can usually be achieved by arranging the task sequence manually, or, alternatively, it can be done using an automated procedure.
in this work we modify the cutting plane mechanism used in the high-dimension function optimization problem to develop an adaptive strategy to automatically perform incremental evolution.
in network inference, solving an easier version of the original task means evolving partial solutions  that can provide some guidance for the search and move toward the target solution gradually.
in the proposed method, parameter sensitivities are calculated and used to determine the priorities of the parameters to be evolved .
the sensitive parameters are selected iteratively, by the proposed sa method and the evolutionary method is employed to search the parameter dimensions accordingly.
an exploration phase is also included in the algorithm for the injection of random effects.
this maintains population diversity and avoids a situation in which individuals move near locally optimal solutions.
the sa and exploration procedures are both performed periodically with specific generation intervals.
in this way, the overall solutions can be derived gradually, and the inferred networks are more robust to internal variations.
the flow of the proposed sa-based incremental evolution approach is illustrated in figure  <dig>  and the details of the approach are described below.
sensitivity-based incremental evolution algorithm()
step 1: initialize the population and start the evolution.
note that steps  <dig> and  <dig> are performed only at certain generations .
step 2: perform sensitivity analysis  to obtain a list of parameters ranked by their sensitivity.
this step is performed only if the evolutionary cycle reaches a pre-defined generation number .
use a threshold  to select the most sensitive parameters.
then update the threshold in order to consider more parameters in the next sensitivity analysis iteration.
set new value ranges for all network parameters.
sensitive parameters are given tight bounds, and insensitive parameters are given loose bounds.
step 3: start the exploration procedure, again only if the evolutionary cycle reaches a pre-defined generation number .
in this phase, a random value is generated for each parameter from its new bounds, replacing each parameter's current value.
step 4: continue the evolutionary computation, using the bounds obtained in step  <dig> to restrict the corresponding parameters.
in the above algorithm, the most critical feature for enabling incremental evolution is found in step  <dig>  the strategy is to perform the sensitivity analysis for network parameters and to specify constraints on them.
a network parameter  is added to the sensitive list and given a pair of tight constraints  if its sensitivity value is less than  the threshold; otherwise, a parameter is added to the insensitive list and given loose constraints.
in the above procedure, the threshold used in step  <dig> to identify the sensitive parameters has an initial value of  <dig>  , and this value is increased gradually with a step value of  <dig> .
the upper and lower bounds of a sensitive parameter in step  <dig> are determined from the best value of this parameter available so far by adding and subtracting an amount equal to twice the value of a small constant .
bounds for the insensitive parameters are defined in a similar way, except that a value of five times the same constant is used.
the dimensions in the sensitive list have higher priorities to be searched, and the corresponding parameter values need to be determined at an earlier evolutionary stage.
as the algorithm describes, the search is accomplished by setting small-range bounds for each sensitive parameter in accordance with the best value found so far, and the evolutionary procedure is performed next to infer a suitable parameter value within the specified bounds.
the creation of constraints for parameters is important in the search for robust solutions, especially when a non-deterministic search method  is employed to derive parameter solutions .
if the parameter value exceeds the specified boundary value, the system dynamics will be influenced, and consequently, the network behaviors will change.
to generate tight constraints for sensitive parameters in the search process is in fact to restrict search regions for the integrated algorithm, which can thus obtain robust networks with desired behaviors.
as indicated in  <cit> , by observing how parameters with varied boundaries direct a system to move toward different system dynamics, we can obtain new insights into a complex biological system and understand the principle of biological design.
as the evolutionary method proceeds, parameters in the insensitive list will be moved to the sensitive list incrementally as the threshold is gradually increased.
therefore, more parameters will be considered sensitive and their values will be determined.
eventually, all the parameter values will be obtained.
because sensitive parameters are more influential to network behavior, enforcing tight constraints on these parameters maintains the robustness of the entire genetic network.
at the same time, it allows other parameters to evolve during the search process.
in this way, the proposed method can infer robust network models with the desired system behavior.
RESULTS
to verify the proposed approach, we conducted a series of experiments on several datasets collected from the literature.
as mentioned, our main goal was to explore the use of sensitivity analysis to infer robust networks, not to compare different computational methods.
therefore, without losing generality, we chose to implement three computational methods, including the traditional pso method, a ga-pso hybrid method, and a differential evolution  method  <cit> , and to build the proposed m-mpsa on them.
the main reason to select the above three methods is that they are the most representative and popular evolutionary algorithms used for optimization tasks.
the same sa strategy can also be embedded into other methods if their implementation details are available.
the first and second phases of the experiments involved comparing the external behavior and internal robustness of the networks inferred by different methods.
four datasets were used.
for each dataset, the above three algorithms were performed with three different settings:  using only the original algorithms;  using our m-mpsa to select sensitive parameters and to derive parameter constraints in the incremental evolution process ; and  using m-mpsa to select sensitive parameters, but with random value bounds during the incremental evolution process.
the third setting was mainly included to demonstrate the importance of using appropriate bounds in the evolutionary process.
in addition, one real dataset was used in the third phase, which made it possible to investigate how the proposed approach can be applied to the study of real gene networks.
this dataset came from a study of gene expression in the sos dna repair system of e. coli.
for this dataset, network models were inferred and the critical parameters were analyzed and discussed.
performance of sensitivity analysis in network modeling
in this set of experiments, four datasets  were used to evaluate the proposed approach.
the first dataset was a five-node artificial network that has become a popular model for comparing different methods in recent studies on network inference  <cit> .
to collect time series data, we started network operations and continued operations for thirty simulation time steps.
the nodes had the following relationships:
Ẋ1= <dig> x3x5- <dig> - <dig> x <dig> Ẋ2= <dig> x <dig> - <dig> x <dig> Ẋ3= <dig> x2- <dig> - <dig> x2- <dig> x <dig>  Ẋ4= <dig> x <dig> x5- <dig> - <dig> x <dig> Ẋ5= <dig> x <dig> - <dig> x <dig>
the second dataset was taken from  <cit> .
it was an eight-gene network created manually by the popular grn simulation software tool genexp.
the third dataset was a ten-node network previously described in  <cit> , given by the following equations:
Ẋ1= <dig> x4x6- <dig> - <dig> x <dig> Ẋ2= <dig> x3x <dig> - <dig> x <dig> Ẋ3= <dig> x1- <dig> x4- <dig> - <dig> x <dig> Ẋ4= <dig> x <dig> x9- <dig> x <dig> Ẋ5= <dig> x <dig> x6- <dig> - <dig> x <dig> Ẋ6= <dig> x <dig> x10- <dig> - <dig> x <dig>  Ẋ7= <dig> x6x10- <dig> - <dig> x <dig> Ẋ8= <dig> x1x2- <dig> x7- <dig> x <dig> Ẋ9= <dig> x3x8- <dig> - <dig> x <dig> Ẋ10= <dig> x <dig> x7- <dig> - <dig> x <dig>
the fourth dataset was part of a real experimental dataset from a study of the rat central nervous system   <cit> .
the original dataset included expression data for  <dig> genes collected at  <dig> time points .
in this experiment, we selected a group of eight nodes representing the largest sub-network of the representative cluster reported in  <cit>  as the target network.
for the above four datasets, three inference algorithms  with three different settings were arranged, and twenty independent runs were conducted for each arrangement.
the population sizes for the four datasets were  <dig>   <dig>   <dig>  and  <dig>  respectively.
tables  <dig>   <dig>   <dig>   <dig> show the results, in which the mean, standard deviation, and best and worst performance of the runs are listed for each arrangement.
fitness results obtained by three inference algorithms with different settings for dataset  <dig>
fitness results obtained by three inference algorithms with different settings for dataset  <dig>
fitness results obtained by three inference algorithms with different settings for dataset  <dig>
fitness results obtained by three inference algorithms with different settings for dataset  <dig>
as shown in tables  <dig>   <dig>   <dig>   <dig>  sa  consistently outperformed the other two settings for inferring gene networks.
it produced the best results for the average, best, and worst fitness values when used in conjunction with any of the three inference methods.
compared with other methods, the ga-pso method produces smaller standard deviations, which indicates that it is more stable than the other methods.
it should also be noted that the performance when using a m-mpsa method with random bounds  was not as good as when using m-mpsa to select parameters and derive bounds for them.
these results confirm the effectiveness of the proposed bound-setting strategy.
evaluation of network robustness
to evaluate the robustness of the networks inferred by three different settings for each algorithm, we compared fitness and the sensitivity values  of the best solutions recorded from the final generations in all runs.
tables  <dig>   <dig>   <dig> list the results of using different settings for the three algorithms.
the values recorded in the tables are the averaged results of the twenty runs performed for each setting.
we can see that results for the four datasets presented in tables  <dig>   <dig>   <dig> are consistent for all three inference algorithms.
they confirm that the proposed sa-based approach  is able to infer networks with lower fitness/error values  and lower sensitivity values  simultaneously.
it should be noted that, because the three algorithms have different computational features , the thresholds  used for constructing the sensitive parameter list in each algorithm were not the same during the evolutionary process.
therefore, the sensitivities obtained from different algorithms are not directly comparable.
comparisons of fitness and sensitivity values of the best solutions collected from the final generations.
experiments were conducted using the pso method.
comparisons of fitness and sensitivity values of the best solutions collected from the final generations.
experiments were conducted using the ga-pso method.
comparisons of fitness and sensitivity values of the best solutions collected from the final generations.
experiments were conducted using the de method.
following the comparisons of the inferred networks, we used an example to further illustrate the effectiveness and importance of using sensitivity analysis to derive appropriate constraints for network parameters in the inference process.
here, we chose the third dataset  and investigated the corresponding parameter correlations.
in this example of eight-gene network, there were  <dig> parameters  to be determined, of which the three most sensitive parameters, p <dig>  p <dig>  and p <dig> , were analyzed.
because the ga-pso algorithm can give the best performance , data were collected from the runs using this method without  and with  the use of the sa method to derive value ranges for network parameters.
the values of the above three parameters recorded from the runs are shown in figure  <dig>  in which the upper and lower parts are the results of the twenty runs conducted for settings  <dig> and  <dig>  respectively.
in this figure, the x-axis indicates the identity of a run, and the y-axis indicates the parameter values.
in the results, we noticed two parameter correlations  in the evolutionary process for each of the experimental settings.
table  <dig> lists the patterns most often obtained from the runs for the two settings .
for setting  <dig>  the two patterns  describe two qualitative relationships: p <dig> >p <dig> >p <dig> and p <dig> >p <dig> >p <dig>  for setting  <dig>  the patterns are p <dig> >p <dig> >p <dig>  and p <dig> >p <dig> >p <dig> .
the runs in which the patterns were observed are also indicated in the tables.
the parameter patterns obtained using two settings.
in table  <dig>  we can see that the qualitative relationship p <dig> >p <dig> >p <dig> appeared in the results for both settings .
after further investigation, we found that this relationship among the three parameters was important for producing good models .
using our specially designed sa method and parameter constraints, this relationship can always be derived and kept in the evolutionary process , but it is not often obtained in the runs with setting  <dig> .
similarly, another parameter relationship, p <dig> >p <dig> >p <dig> , obtained from setting  <dig>  was useful for inferring good models.
to observe how the three parameters varied during the runs, see figures  <dig>   <dig>  these figures illustrate typical examples of the patterns presented in table  <dig>  in the figures, the x-axis indicates the generation number , and the y-axis indicates the parameter value.
as can be clearly seen in figures  <dig>   <dig> the parameter values changed actively to move toward the appropriate positions to obtain a successful model in the runs using setting  <dig>  whereas the values remained static in the runs using setting  <dig>  these results verify the effectiveness of using sa in network inference.
evaluation of the proposed approach on a real dataset
after evaluating the performance of the proposed sa-based approach in network inference, we conducted a set of experiments to investigate how our approach can be applied to the study of a real gene network.
because the ga-pso algorithm has been shown to outperform other methods in the above experiments, we used this method with two settings  to conduct experiments on a real dataset.
the data set used in this experimental phase comes from a study of the sos dna repair system in e. coli.
figure  <dig> illustrates the gene regulation that occurs in this system  <cit> .
in this system, the lexa protein  maintains its expression level in a normal state and is bound to the promoter regions of sos genes, including uvrd, umud, lexa, uvra, reca, and polb.
when dna damage occurs, the reca protein senses the damage and mediates lexa autocleavage.
the decrease in lexa relieves the repression of the sos genes.
the expression of these genes then activates the sos repair system.
once the damage has been repaired, the concentration of reca drops and this gene stops mediating lexa autocleavage.
the lexa level increases and begins to repress the sos genes again.
inference of expression profiles with the decomposed s-system model
in the experiments, six genes were selected from the original experimental data reported in  <cit> .
the genes selected were uvrd, umud, lexa, uvra, reca, and polb.
these were selected because the interactions of the 6-gene network have been well studied and commonly used in related studies, and the corresponding network has been inferred successfully  <cit> .
using the regulatory relationships that have been described, we can validate our proposed network inference method by comparing our results to those reported previously.
though there have been several studies on inferring the sos repair system, none of them used a tightly coupled s-system model  to represent the gene network, and none of them inferred such a model from expression data.
to the best of our knowledge, the most relevant work is  <cit> , in which the authors utilized a decoupled s-system model to infer the sos repair system.
in a decoupled s-system, a tightly coupled system of non-linear differential equations is decomposed and analyzed as several differential equations  <cit> , each of which can describe a specific gene and can then be separately inferred .
motivated by this research, we thus adopted decoupled differentials to describe gene profiles, and we employed the proposed approach to infer gene interactions.
in this way, we can not only examine the inferred network behaviors, but also compare the inferred gene regulations to the ones obtained in the other studies mentioned above.
the sos repair system dataset included  <dig> sampling points for each gene.
the six major genes  were inferred through the decoupled s-system model.
in the experiments, twenty independent runs were conducted for each of the two settings.
table  <dig> presents the results for the six genes of the sos repair system dataset obtained by using settings  <dig> and  <dig>  the mean, variance, best and worst fitness, and average parameter sensitivity values of the twenty runs are listed.
again, the results show that the sa-based approach  outperformed the original algorithm in terms of both fitting the external behavior  and exploring the internal robustness  of a gene network.
in addition, to observe the inferred network behaviors, figures  <dig>   <dig> compare the inferred and target behaviors.
in these figures, the x-axis represents time, and the y-axis represents the concentrations of particular gene products.
as shown in the figures, network behaviors very similar to those of the real system can be inferred by our approach.
results obtained by the ga-pso algorithm with two settings.
analysis of critical parameters of the sos repair system
to investigate the critical parameters that have significant influences on the system dynamics of the sos repair system, we summarized the most sensitive parameters identified by the proposed approach in table  <dig>  there were  <dig> parameters selected and marked as crucial interactions, each of which represented a regulatory relationship between two genes .
the letters g and h in the table indicate which of the kinetic orders gi,j and hi,j of gene i  is selected.
as mentioned previously, the identified parameters  can also be used to determine the network structure of the system to be inferred.
by exploiting the structural information, the evolutionary algorithm  can infer better models with robust system dynamics during the search process.
in the sos case, the  <dig> parameters indicated in table  <dig> can be used to derive the scaffold of the sos network, as shown in figure  <dig>
the most sensitive parameters identified by the proposed approach.
the letter g indicates that gi,j is the most influential parameter for gene i; while the letter h indicates that hi,j is the most influential parameter.
the asterisk means the regulation known in the literature .
to compare the gene regulation events identified by our approach to those found in the literature, we highlighted the gene interactions in table  <dig>  and summarized the results in table  <dig>  the results show that, of the thirteen gene regulation events found in our experiments, ten matched regulatory relationships known from other studies.
for example, in figure  <dig>  the negative regulation of lexa, uvra, and uvrd by lexa has been successfully identified as the most crucial interaction in determining system dynamics.
still, it is notable that the known negative regulation of lexa by reca was not recognized.
the main reason for this is that our approach chose sensitive parameters from the kinetic orders of the s-system model on the gene level, but the above relationship was in fact interpreted as a regulation of protein lexa by protein reca  <cit> .
the concentration of protein reca depends on the sensing of dna damage, and the interaction between reca and lexa therefore depends on the events occurring between dna damage and dna repair.
because the on or off state of dna damage sensing was not directly described in the s-system model, the regulation of lexa by reca on the protein level was not observable on the basis of kinetic orders.
in other words, the absence of this regulation was due to the model representation, rather than the inference approach.
this shows that, to infer the sos repair system accurately on both the protein and gene levels, a more comprehensive mathematical model is needed to describe the regulatory details.
with only one exception  that was caused by limitations of the model itself, our approach has shown its strength in inferring real gene system that show the expected network behavior and internal network robustness.
summary of critical gene regulatory relationships obtained by the proposed framework with respect to the results in the literature.
the parameters gi,j and hi,j are kinetic orders of the s-system.
after identifying the most sensitive parameters of the sos repair system, we again analyzed the gene regulation relationships obtained from the runs to observe the correlation of parameters.
the importance of investigating the relationships among parameter values of genes has been emphasized in the study of systems biology, especially when the boundary of a parameter value is taken into consideration  <cit> .
as presented and discussed in the second experimental phase, the qualitative relationships among genes can be regarded as special patterns  that determined the system dynamics of a candidate model.
by exploiting these patterns, the inference algorithm can refine the values of kinetic orders to obtain better solutions.
in the case of the sos repair system, the patterns of regulatory interactions can be derived as in other cases, by recording and analyzing the values of the most sensitive parameters.
taking gene lexa as an example , the algorithm categorized genes lexa and uvra as the critical gene regulators for lexa and two regulatory relationships, lexa -| lexa and uvra -| lexa can be established on the basis of this information.
figure  <dig> illustrates the parameter values related to the patterns derived from the twenty runs.
in this figure, the x-axis indicates the run number, and the y-axis represents the parameter values of the specific pattern  obtained in each run.
in this figure, the upper and lower parts present the results of inferring the network without  and with  using the sa method, respectively.
the regulation patterns for gene lexa are summarized in table  <dig>  as shown in table  <dig>  for experimental runs using setting  <dig>  fifteen out of twenty runs consistently produced the same pattern uvra -| lexa > lexa -| lexa .
this pattern can guide the inference algorithm to find better solutions.
in contrast, the qualitative relationships between the two relationships uvra -| lexa and lexa -| lexa obtained from the runs using setting  <dig> are inconsistent and unstable .
the above results show that the proposed approach can derive parameter bounds  to restrict the variation of lexa, and it can therefore infer models with better fitness and more stable  system dynamics.
without using the sa method, the inference algorithm alone cannot guarantee the preservation of useful patterns during the iterative inference process.
pattern-ns1:
pattern-ns2:
pattern-s1:
the patterns of gene regulations found for target gene lexa.
tables  <dig>   <dig>   <dig> list the regulation patterns analyzed for four other target genes in the sos system: lexa, uvra, uvrd, and umud.
figures  <dig>   <dig>   <dig> show the parameter values involved.
the genes reca and polb were ignored because there was only one gene regulatory relationship found for each of them .
this means that no qualitative relationship can be established for these genes.
the results in tables  <dig>   <dig>   <dig> are consistent with table  <dig>  which again shows the success of using sa in network inference.
pattern-ns1:
pattern-s1:
the patterns of gene regulations found for target gene uvra.
pattern-ns1:
pattern-s1:
pattern-s2:
the patterns of gene regulations found for target gene uvrd.
pattern-ns1:
pattern-ns2:
pattern-s1:
pattern-s2:
the patterns of gene regulations found for target gene umud.
CONCLUSIONS
in this work, we developed a sensitivity-based incremental evolution approach to cope with one important issue, network robustness, which has not been addressed in gene network inference.
our approach included two parts.
the first part was a sensitivity analysis method that was used to select sensitive parameters for deriving value bounds of these parameters.
the second part was an evolutionary algorithm that took the bounds as constraints to perform parameter optimization.
this process leads to inferred networks that are robust and produce the desired behaviors.
to validate the proposed approach, a series of experiments were conducted to evaluate the external behaviors and internal robustness of the networks inferred by different methods.
the results show that the proposed sa-based approach outperformed other methods in all datasets.
moreover, we analyzed in detail the results obtained from real time-series expression data for the sos repair system.
the analyses indicate that our approach can identify the critical parameters and use them to establish regulatory relationships among genes.
by enforcing these relationships in the repetitive search process, our approach can successfully infer robust networks.
competing interests
the authors declare that they have no competing interests.
authors' contributions
yh undertook the experimental implementation, made statistical analysis and wrote a part of the manuscript.
wl conceived the project, designed the algorithm, wrote a part of the manuscript, and made modifications as well as final revisions.

 gene regulatory networks are essential for controlling cellular metabolism and the organismal development. under the command of transcription factors, each gene influences the activity of the cell by generating messenger rna that guides the synthesis of proteins by ribosomes in the cytoplasm. gene network modeling uses gene expression data to characterize the phenotypic behavior of a system under study. with the reconstructed networks, biologists can generate and test hypotheses to further understand the complex phenomena that occur in nature systems and to explore the dynamics of those systems. modeling grns manually on the basis of the experimentally measured time-series data takes a considerable amount of time. therefore, an automated reverse-engineering procedure is recommended. this procedure involves altering the gene network in some way, observing the outcome, and using computational methods to infer the underlying principles of the network. to derive a realistic model, available domain knowledge can be integrated into the computational methods. figure illustrates the general procedure of inferring grns from quantitative expression data. the inferring/modeling block of this figure indicates the computational procedure used to derive network parameters for a given model, to build and simulate the model, and to evaluate it by comparing the behavior of the inferred model with the original data set. in addition to making use of expression data, one recently developed strategy combines information from various sources to narrow down the search space in the network. this strategy shortens the time and effort required for the validation and discovery of networks. for example, if some gene names are known then they can be mapped into knowledge bases to extract biological knowledge that can be used to determine the network structures. our goal in this study is to establish a methodology for network inference and to investigate aspects of network dynamics that have not been addressed previously. to infer a network with desired system behavior, the most important steps are to select a network model and then to fit the network's structural parameters to the available expression data. many models have been proposed to address different levels of biological details, ranging from the very abstract to the very concrete. abstract models are easy to simulate and therefore less computationally taxing, but it has been proven that they are not able to capture certain system behaviors. in contrast, concrete models are more suitable for simulating biochemical process realistically, but, because of their computational complexities, these models can only be applied to small systems. to capture the underlying physical phenomena of a gene network we took one of the most popular and well-studied concrete models, the s-system model, as an example of a biological network, and we used both simulated and collected real gene expression data to reconstruct the model. the s-system model is a type of ordinary differential equation model. it consists of a particular set of tightly coupled odes in which the component processes are characterized by power law functions. the system structure of an s-system is described by the following equation: dxidt=αi ∏j=1nxjgi,j-βi ∏j=1nxjhi,j here, xi is the expression level of gene i, and n is the number of genes in a genetic network. the non-negative parameters αi and βi are rate constants that indicate the direction of mass flow. the real-number exponents gi,j and hi,j are kinetic orders that reflect the strength of interactions from genes j to i. the above set of parameters defines an s-system model. to infer an s-system model is, thus, to estimate all of the 2n parameters simultaneously. although non-linear odes can more accurately model the system dynamics of gene networks, they are difficult to solve by traditional optimization techniques. to reduce the complexity of such models, di bernardo et. al developed an approach that involves taking a series of steady-state gene expression measurements following transcriptional perturbations and reducing the linear odes to a linear regression problem that can be solved with relevant techniques. in addition, some researchers employed the strategies of gene clustering for dimension reduction and gene classification for identification of expression patterns to reduce task complexity. as evolutionary algorithms have been widely used to solve many difficult optimization tasks with good results, they have been suggested as a way to infer gene networks. among many studies on this topic, the most relevant are those involving the use of eas to infer s-system models. in these studies, the network parameters of the s-system were arranged as a string of floating numbers that could be evolved by genetic operators. in recent studies of ea-based parameter inference, one critical problem, network robustness, has not been addressed. it is important to investigate the effect of network parameter perturbations on the overall system. these parameters not only are the numerical components of the model, but also represent the activities or interactions among proteins, transcription factors, and mrnas in the gene network. each network parameter has an important role in determining the behavior of a biological system. parameters that are very sensitive to variation can introduce fragility into the system. in fact, recent research has found that parameters can define the dynamics of a model, and more important, studies have shown that enforcing constraints on the parameters can limit model dynamics. specifically, by measuring and analyzing the variations of network parameters and their effects on relevant genes, the gene-gene interactions of an inferred model can be interpreted. to ensure the robustness of the inferred network and to further investigate gene interactions, it is important to derive an acceptable value range for each parameter and to restrict the parameter's value to the specified range during the network reconstruction process. therefore, we take parameter sensitivity into consideration in the network reconstruction procedure so that robust results can be obtained. sensitivity analysis is an indispensable technique that can help researchers to investigate the parameter robustness properties of an inferred network. by varying the parameter values within a certain range and performing statistical calculations to measure the system instability, researchers can identify the critical parameters or discover seemingly unimportant parameters that may have a positive or negative influence on a network. sa techniques have been used in many biological studies, such as genetic circuit design, mammalian circadian clock modeling, and target prediction in signaling pathways. in general, the sensitivity of a parameter is defined as: spm=∂m/m∂p/p=percentage change in mpercentage change in p here, p represents the parameters that are varied in a given range, m is the mathematical function describing the system behavior, and ∂m means the change in m due to the value changed in ∂p with respect to p. in a dynamic system, there are two types of methods often used to conduct sa: local sa and global sa. local sa can measure the sensitivity of a parameter from a given range even if the system structure is unknown, but it only considers one parameter at a time and ignores the interactions between parameters. alternatively, global sa can examine parameter interactions with different parameter magnitudes simultaneously, but this can only be done if the system structure is known in advance. the above two methods have been compared extensively. it should be noted that neither local sa nor global sa is superior: the choice of method depends on the specific application and the information available. in the biological cases that involve the use of sa techniques, the sensitivity of parameters should be investigated from a biological point of view in order to validate the results. from the perspective of parameter dynamics, we propose a new approach to infer robust networks. our approach includes two major parts. the first part is a sensitivity analysis procedure that selects sensitive network parameters, determines value ranges for them, and then sends the parameters with their constraints to the inference mechanism. the second part is an evolutionary method responsible for inferring networks from the results obtained in the first part. this approach, to our knowledge, is the first work in network modeling that integrates sensitivity analysis into an inference algorithm to consider both internal and external characteristics of the inferred network. to validate the proposed approach, a series of experiments has been conducted, in which the proposed sa procedure was coupled with different eas. the results have been analyzed, and they show that this approach can be used to infer robust networks successfully from artificial and real gene expression profiles. methods parameter sensitivity analysis in the process of inferring a gene network from expression data the genes interact with each other, and the network structure is generally unclear during the modeling process. most studies of global sa use known pathways as examples and can therefore easily select the critical parameters and control the total number of parameters for further analysis. in the type of network inference problem addressed here, it is not possible to identify and select the most important parameters in order to perform global sa, because no prior information about network structure is available. to consider multiple network parameters simultaneously in such situations, we present a new method that is a modification of a widely used global sa technique, called multi-parameter sensitivity analysis. mpsa is based on monte carlo simulations and on quantitative comparisons of cumulative frequency via corresponding pearson correlation coefficients among parameters to identify those that are most sensitive. it can be used to explore critical genetic reactions and to examine the robustness of a gene network as other global sa methods do. in the method described here, the objective function of the sa method is the fitness function that is to be optimized in the evolutionary algorithm, which is the mean squared error over the time course: ∑i=1n∑t=1txia-xidxid in this expression, xid is the desired expression level of gene i at time t, xia is the actual value obtained from the inferred model, n is the number of genes in the network, and t is the number of time points used to measure gene expression data during the period. our modified method, m-mpsa, inherits qualities from mpsa but has several advantages. first, mpsa is easy to implement and modify. it allows users to select one parameter at a time or to consider multiple parameters at the same time to perform sensitivity analysis. this characteristic is crucial because, when the network structure is unknown or uncertain, selecting one parameter at a time is a good way to reduce the complexity of the search space. researchers can then concentrate on each parameter and temporarily neglect the combination sets at the selection phase. second, compared with most of the global sa methods, the mathematical equations of mpsa are more manifest and succinct, and the parameter sensitivity is thus easy to calculate. finally, the value ranges of parameters used in mpsa are the same as in an ode model. this means that the sensitivity of each parameter can be directly mapped into an ode model without any extra computation. in previous studies, researchers have shown that the most influential parameters play key roles in modeling a grn; that is, a model's behavior varies largely because of its sensitive parameters. therefore, our current study draws on the characteristics of sensitive parameters, and it incorporates the sensitivity analysis method into the incremental evolution approach to infer a grn. the proposed m-mpsa method is described in the sensitivity analysis algorithm below. it uses an iterative process to calculate the sensitivity of each parameter, and it then ranks the sensitivities of all parameters. the input of this algorithm is the set of network parameters to be determined; and the output, a list of parameters ranked by their sensitivity values. with the parameter sensitivities ranked, the evolutionary algorithm can then infer robust solutions by exploiting the sensitive parameters that are likely to significantly influence the genetic model as a whole. sensitivity analysis algorithm() step 1: select one parameter, i, at a time.. step 2: set the parameter range ri for each s-system parameter, i. initially, the commonly used search regions and are taken as the initial settings for the parameters of kinetic orders and rate constants respectively, as suggested in. these settings are used to define the parameter range ri, by adding/subtracting an amount of one third of the initial settings to/from the current value of parameter i. step 3: for each ri, a set of random values are created with a uniform distribution within the specified range. each random value for parameter i, together with the other parameters, constitutes a candidate solution and calculate the objective function value for each of the random values generated in step step 5: determine whether the objective function value of each random point obtained in step is acceptable or unacceptable by a given threshold cr. the parameter value of a random point is classified as unacceptable if its corresponding objective function value is greater than cr; otherwise, the parameter value is acceptable. go back to step until all parameters have been dealt with. step 6: calculate the sensitivity for parameter i by using its cumulative frequency. cfi measures the similarity of two statistical distributions that are formed by the acceptable and unacceptable values produced from steps 3- for parameter i. the parameters are ranked according to their cf, and the ones with relatively low cf are considered to be sensitive. step 7: output a parameter list in which the parameters are sorted by their cf values. in addition to presenting the m-mpsa algorithm, we provide here a walkthrough example to illustrate how it operates in practice. in the example, we randomly choose two parameters from one of the datasets used in the experimental section, and we calculate the sensitivity values of these parameters. the dataset consisted of ten nodes, and there were altogether parameters to be estimated. without losing generality, we assume that the algorithm selected the parameter p, and then set the range r for p. as mentioned above, we used one third of the original search region to set up the bounds for each parameter. hence, if p had a value of, then the lower bound for this parameter was - )), and the upper bound was. then, random values were generated within the range for p, and their corresponding objective function values were calculated. next, a threshold cr was determined to categorize the above evaluation results as acceptable or unacceptable. if the current best result was, then a threshold of was used to determine whether a result was acceptable or not. after that, two statistical distributions were constructed for the acceptable and unacceptable results identified above. for statistical purposes, the value range was divided into ten intervals, and for each interval the numbers of acceptable and unacceptable values were counted. figure shows the distributions for p in which acceptable and unacceptable numbers were marked. the cumulative frequency distributions were then constructed accordingly, as presented in figure. finally, the sensitivity value of p was obtained by measuring the correlation coefficients of the two cumulative frequency distributions. in this case, the sensitivity value of p was. a similar procedure was performed for the other parameter, p and a sensitivity value of was obtained. because p had a sensitivity value lower than that of p it was considered to be more sensitive than p in this way, after sensitivity values of all parameters were calculated, they were ranked and sent to the incremental evolution algorithm. as mentioned above, more sensitive parameters are more likely to have a significant influence on the entire genetic model. therefore, if no prior knowledge can be applied to the reconstruction of a large grn, researchers can conduct sensitivity analysis for parameter selection to guide the network inference method, and observe how numerical quantities within different intervals can formulate candidate results. taking p as an example, the acceptable values for this parameter fell in the range that consisted of acceptable network models in total. we can exploit this result to infer p and explore the relationships of this parameter with other genetic reactions. in fact, researchers have shown that the most influential parameters play important roles in modeling a grn. that is, the behavior variation of a network model largely depends on its sensitive parameters. in figure we can see that the variation of the parameter value range influenced the system dynamics. that is, the system robustness was significantly affected by the value range specified for a parameter. for instance, the two distributions for p in figure show that, if the value range of p changed from to, then the acceptable values for p changed from to this means that the inference algorithm can maintain the system robustness if it can locate p within the interval. the above analysis concludes that the proposed sa method can be used to find sensitive parameters and to derive suitable value ranges as search constraints, which can then be used in the inference algorithm to construct robust networks. evolutionary algorithm for parameter optimization to derive the network parameters of s-system model, we also implement an evolutionary algorithm for parameter optimization to work with the above sa method. ea is a population-based approach that evaluates many solutions simultaneously in the search space, and it is likely to find a global solution for a given problem. recently, a new population-based optimization technique, particle swarm optimization, was proposed as an alternative to the traditional eas. pso tries to mimic the goal-seeking behavior of biological swarms. the standard pso algorithm contains a set of particles and operates in an iterative manner. each particle is characterized by its position and velocity, and it moves in the solution space. the position of each particle represents a potential solution that is evaluated by a predefined fitness function. pso has some attractive characteristics. in particular, it has memory, so that knowledge of a good solution can be retained by all particles. during the iterative search process, each particle remembers its previous best position and the best position of any particle in the swarm. then the particle uses the position information to modify its position and velocity, and continues its movement in the search space. the details of the pso algorithm can be found in. some performance comparisons between pso and the most popular ea method, the genetic algorithm, have been made, underscoring the reliability and the convergence speed of both methods. however, the result tends to be inconclusive: each of the algorithms has shown better performance than the other for some particular applications. consequently, hybrid techniques were proposed to effectively exploit the qualities and the uniqueness of the two methods. it is now commonly agreed that the hybrid methods can lead to further performance improvements. therefore, in this work, we develop a hybrid method to exploit the solution memory of pso and the genetic operations of ga for network inference. as has been mentioned, our primary goal is to investigate how parameter sensitivity can be used to guide an ea and derive robust networks. to concentrate on this issue, we have chosen to implement a popular ga-pso method for parameter optimization, and we have built our sa method upon it. breeding swarms is a frequently used ga-pso hybrid method that has demonstrated its good performance in many benchmark test functions. therefore, we adopt this method for our network inference application. the hybrid ga-pso method used in this work is illustrated in figure initially, a population is randomly generated and evaluated. the individual solutions are ranked according to their fitness values, and then they are divided into two parts for running the pso and ga methods separately. as shown in figure the best individuals % of the whole population) are preserved and enhanced by the pso procedure. they are then sent to the next generation. meanwhile, on the right-hand side of the figure, the remaining individuals are discarded. to replace the removed individuals, a tournament selection scheme is used to choose the same number of individuals from the best ones, and the selected individuals are used to create new individuals by the ga procedure. in our implementation, we take a direct encoding scheme to represent solutions for both the ga and pso parts, in which the network parameters are arranged as a linear string chromosome of floating-point numbers. the goal here is to minimize the accumulated discrepancy between the gene expression data recorded in the dataset and the values produced by the inferred model. therefore, the error function defined previously is used directly for performance measurement. for the pso part, the equations for updating the particle's velocity and position are the same as the ones listed in the original pso work, and for the ga part, the operations of crossover and mutation described in are used. sensitivity-based incremental evolution method for network inference as described in the first section, to infer the s-system model is to determine the 2n parameters simultaneously. solving this high-dimensional problem is difficult, especially when the complexity of regulation increases along with the number of genes involved. one promising approach to this problem is to adopt the concept of incremental evolution to infer the large number of parameters. the underlying principle of incremental evolution is that a population is first evolved to solve an easier version t' of the original complex task t, in which the solution region of t is more accessible from region t'. more task versions with incremental complexity can be arranged so that the original task can be achieved progressively. evidently, the main task involved in implementing incremental evolution is the formulation of a scheme to transfer the goal task into another more evolvable task. in the process of task transformation, the structure of the environment and the overall goal must be preserved. this can usually be achieved by arranging the task sequence manually, or, alternatively, it can be done using an automated procedure. in this work we modify the cutting plane mechanism used in the high-dimension function optimization problem to develop an adaptive strategy to automatically perform incremental evolution. in network inference, solving an easier version of the original task means evolving partial solutions that can provide some guidance for the search and move toward the target solution gradually. in the proposed method, parameter sensitivities are calculated and used to determine the priorities of the parameters to be evolved. the sensitive parameters are selected iteratively, by the proposed sa method and the evolutionary method is employed to search the parameter dimensions accordingly. an exploration phase is also included in the algorithm for the injection of random effects. this maintains population diversity and avoids a situation in which individuals move near locally optimal solutions. the sa and exploration procedures are both performed periodically with specific generation intervals. in this way, the overall solutions can be derived gradually, and the inferred networks are more robust to internal variations. the flow of the proposed sa-based incremental evolution approach is illustrated in figure and the details of the approach are described below. sensitivity-based incremental evolution algorithm() step 1: initialize the population and start the evolution. note that steps and are performed only at certain generations. step 2: perform sensitivity analysis to obtain a list of parameters ranked by their sensitivity. this step is performed only if the evolutionary cycle reaches a pre-defined generation number. use a threshold to select the most sensitive parameters. then update the threshold in order to consider more parameters in the next sensitivity analysis iteration. set new value ranges for all network parameters. sensitive parameters are given tight bounds, and insensitive parameters are given loose bounds. step 3: start the exploration procedure, again only if the evolutionary cycle reaches a pre-defined generation number. in this phase, a random value is generated for each parameter from its new bounds, replacing each parameter's current value. step 4: continue the evolutionary computation, using the bounds obtained in step to restrict the corresponding parameters. in the above algorithm, the most critical feature for enabling incremental evolution is found in step the strategy is to perform the sensitivity analysis for network parameters and to specify constraints on them. a network parameter is added to the sensitive list and given a pair of tight constraints if its sensitivity value is less than the threshold; otherwise, a parameter is added to the insensitive list and given loose constraints. in the above procedure, the threshold used in step to identify the sensitive parameters has an initial value of, and this value is increased gradually with a step value of. the upper and lower bounds of a sensitive parameter in step are determined from the best value of this parameter available so far by adding and subtracting an amount equal to twice the value of a small constant. bounds for the insensitive parameters are defined in a similar way, except that a value of five times the same constant is used. the dimensions in the sensitive list have higher priorities to be searched, and the corresponding parameter values need to be determined at an earlier evolutionary stage. as the algorithm describes, the search is accomplished by setting small-range bounds for each sensitive parameter in accordance with the best value found so far, and the evolutionary procedure is performed next to infer a suitable parameter value within the specified bounds. the creation of constraints for parameters is important in the search for robust solutions, especially when a non-deterministic search method is employed to derive parameter solutions. if the parameter value exceeds the specified boundary value, the system dynamics will be influenced, and consequently, the network behaviors will change. to generate tight constraints for sensitive parameters in the search process is in fact to restrict search regions for the integrated algorithm, which can thus obtain robust networks with desired behaviors. as indicated in, by observing how parameters with varied boundaries direct a system to move toward different system dynamics, we can obtain new insights into a complex biological system and understand the principle of biological design. as the evolutionary method proceeds, parameters in the insensitive list will be moved to the sensitive list incrementally as the threshold is gradually increased. therefore, more parameters will be considered sensitive and their values will be determined. eventually, all the parameter values will be obtained. because sensitive parameters are more influential to network behavior, enforcing tight constraints on these parameters maintains the robustness of the entire genetic network. at the same time, it allows other parameters to evolve during the search process. in this way, the proposed method can infer robust network models with the desired system behavior.  to verify the proposed approach, we conducted a series of experiments on several datasets collected from the literature. as mentioned, our main goal was to explore the use of sensitivity analysis to infer robust networks, not to compare different computational methods. therefore, without losing generality, we chose to implement three computational methods, including the traditional pso method, a ga-pso hybrid method, and a differential evolution method, and to build the proposed m-mpsa on them. the main reason to select the above three methods is that they are the most representative and popular evolutionary algorithms used for optimization tasks. the same sa strategy can also be embedded into other methods if their implementation details are available. the first and second phases of the experiments involved comparing the external behavior and internal robustness of the networks inferred by different methods. four datasets were used. for each dataset, the above three algorithms were performed with three different settings: using only the original algorithms; using our m-mpsa to select sensitive parameters and to derive parameter constraints in the incremental evolution process; and using m-mpsa to select sensitive parameters, but with random value bounds during the incremental evolution process. the third setting was mainly included to demonstrate the importance of using appropriate bounds in the evolutionary process. in addition, one real dataset was used in the third phase, which made it possible to investigate how the proposed approach can be applied to the study of real gene networks. this dataset came from a study of gene expression in the sos dna repair system of e. coli. for this dataset, network models were inferred and the critical parameters were analyzed and discussed. performance of sensitivity analysis in network modeling in this set of experiments, four datasets were used to evaluate the proposed approach. the first dataset was a five-node artificial network that has become a popular model for comparing different methods in recent studies on network inference. to collect time series data, we started network operations and continued operations for thirty simulation time steps. the nodes had the following relationships: Ẋ1= x3x5- - x Ẋ2= x - x Ẋ3= x2- - x2- x  Ẋ4= x x5- - x Ẋ5= x - x the second dataset was taken from. it was an eight-gene network created manually by the popular grn simulation software tool genexp. the third dataset was a ten-node network previously described in, given by the following equations: Ẋ1= x4x6- - x Ẋ2= x3x - x Ẋ3= x1- x4- - x Ẋ4= x x9- x Ẋ5= x x6- - x Ẋ6= x x10- - x  Ẋ7= x6x10- - x Ẋ8= x1x2- x7- x Ẋ9= x3x8- - x Ẋ10= x x7- - x the fourth dataset was part of a real experimental dataset from a study of the rat central nervous system. the original dataset included expression data for genes collected at time points. in this experiment, we selected a group of eight nodes representing the largest sub-network of the representative cluster reported in as the target network. for the above four datasets, three inference algorithms with three different settings were arranged, and twenty independent runs were conducted for each arrangement. the population sizes for the four datasets were and respectively. tables show the results, in which the mean, standard deviation, and best and worst performance of the runs are listed for each arrangement. fitness results obtained by three inference algorithms with different settings for dataset fitness results obtained by three inference algorithms with different settings for dataset fitness results obtained by three inference algorithms with different settings for dataset fitness results obtained by three inference algorithms with different settings for dataset as shown in tables sa consistently outperformed the other two settings for inferring gene networks. it produced the best results for the average, best, and worst fitness values when used in conjunction with any of the three inference methods. compared with other methods, the ga-pso method produces smaller standard deviations, which indicates that it is more stable than the other methods. it should also be noted that the performance when using a m-mpsa method with random bounds was not as good as when using m-mpsa to select parameters and derive bounds for them. these results confirm the effectiveness of the proposed bound-setting strategy. evaluation of network robustness to evaluate the robustness of the networks inferred by three different settings for each algorithm, we compared fitness and the sensitivity values of the best solutions recorded from the final generations in all runs. tables list the results of using different settings for the three algorithms. the values recorded in the tables are the averaged results of the twenty runs performed for each setting. we can see that results for the four datasets presented in tables are consistent for all three inference algorithms. they confirm that the proposed sa-based approach is able to infer networks with lower fitness/error values and lower sensitivity values simultaneously. it should be noted that, because the three algorithms have different computational features, the thresholds used for constructing the sensitive parameter list in each algorithm were not the same during the evolutionary process. therefore, the sensitivities obtained from different algorithms are not directly comparable. comparisons of fitness and sensitivity values of the best solutions collected from the final generations. experiments were conducted using the pso method. comparisons of fitness and sensitivity values of the best solutions collected from the final generations. experiments were conducted using the ga-pso method. comparisons of fitness and sensitivity values of the best solutions collected from the final generations. experiments were conducted using the de method. following the comparisons of the inferred networks, we used an example to further illustrate the effectiveness and importance of using sensitivity analysis to derive appropriate constraints for network parameters in the inference process. here, we chose the third dataset and investigated the corresponding parameter correlations. in this example of eight-gene network, there were parameters to be determined, of which the three most sensitive parameters, p p and p, were analyzed. because the ga-pso algorithm can give the best performance, data were collected from the runs using this method without and with the use of the sa method to derive value ranges for network parameters. the values of the above three parameters recorded from the runs are shown in figure in which the upper and lower parts are the results of the twenty runs conducted for settings and respectively. in this figure, the x-axis indicates the identity of a run, and the y-axis indicates the parameter values. in the results, we noticed two parameter correlations in the evolutionary process for each of the experimental settings. table lists the patterns most often obtained from the runs for the two settings. for setting the two patterns describe two qualitative relationships: p >p >p and p >p >p for setting the patterns are p >p >p and p >p >p. the runs in which the patterns were observed are also indicated in the tables. the parameter patterns obtained using two settings. in table we can see that the qualitative relationship p >p >p appeared in the results for both settings. after further investigation, we found that this relationship among the three parameters was important for producing good models. using our specially designed sa method and parameter constraints, this relationship can always be derived and kept in the evolutionary process, but it is not often obtained in the runs with setting. similarly, another parameter relationship, p >p >p, obtained from setting was useful for inferring good models. to observe how the three parameters varied during the runs, see figures these figures illustrate typical examples of the patterns presented in table in the figures, the x-axis indicates the generation number, and the y-axis indicates the parameter value. as can be clearly seen in figures the parameter values changed actively to move toward the appropriate positions to obtain a successful model in the runs using setting whereas the values remained static in the runs using setting these results verify the effectiveness of using sa in network inference. evaluation of the proposed approach on a real dataset after evaluating the performance of the proposed sa-based approach in network inference, we conducted a set of experiments to investigate how our approach can be applied to the study of a real gene network. because the ga-pso algorithm has been shown to outperform other methods in the above experiments, we used this method with two settings to conduct experiments on a real dataset. the data set used in this experimental phase comes from a study of the sos dna repair system in e. coli. figure illustrates the gene regulation that occurs in this system. in this system, the lexa protein maintains its expression level in a normal state and is bound to the promoter regions of sos genes, including uvrd, umud, lexa, uvra, reca, and polb. when dna damage occurs, the reca protein senses the damage and mediates lexa autocleavage. the decrease in lexa relieves the repression of the sos genes. the expression of these genes then activates the sos repair system. once the damage has been repaired, the concentration of reca drops and this gene stops mediating lexa autocleavage. the lexa level increases and begins to repress the sos genes again. inference of expression profiles with the decomposed s-system model in the experiments, six genes were selected from the original experimental data reported in. the genes selected were uvrd, umud, lexa, uvra, reca, and polb. these were selected because the interactions of the 6-gene network have been well studied and commonly used in related studies, and the corresponding network has been inferred successfully. using the regulatory relationships that have been described, we can validate our proposed network inference method by comparing our results to those reported previously. though there have been several studies on inferring the sos repair system, none of them used a tightly coupled s-system model to represent the gene network, and none of them inferred such a model from expression data. to the best of our knowledge, the most relevant work is, in which the authors utilized a decoupled s-system model to infer the sos repair system. in a decoupled s-system, a tightly coupled system of non-linear differential equations is decomposed and analyzed as several differential equations, each of which can describe a specific gene and can then be separately inferred. motivated by this research, we thus adopted decoupled differentials to describe gene profiles, and we employed the proposed approach to infer gene interactions. in this way, we can not only examine the inferred network behaviors, but also compare the inferred gene regulations to the ones obtained in the other studies mentioned above. the sos repair system dataset included sampling points for each gene. the six major genes were inferred through the decoupled s-system model. in the experiments, twenty independent runs were conducted for each of the two settings. table presents the results for the six genes of the sos repair system dataset obtained by using settings and the mean, variance, best and worst fitness, and average parameter sensitivity values of the twenty runs are listed. again, the results show that the sa-based approach outperformed the original algorithm in terms of both fitting the external behavior and exploring the internal robustness of a gene network. in addition, to observe the inferred network behaviors, figures compare the inferred and target behaviors. in these figures, the x-axis represents time, and the y-axis represents the concentrations of particular gene products. as shown in the figures, network behaviors very similar to those of the real system can be inferred by our approach. results obtained by the ga-pso algorithm with two settings. analysis of critical parameters of the sos repair system to investigate the critical parameters that have significant influences on the system dynamics of the sos repair system, we summarized the most sensitive parameters identified by the proposed approach in table there were parameters selected and marked as crucial interactions, each of which represented a regulatory relationship between two genes. the letters g and h in the table indicate which of the kinetic orders gi,j and hi,j of gene i is selected. as mentioned previously, the identified parameters can also be used to determine the network structure of the system to be inferred. by exploiting the structural information, the evolutionary algorithm can infer better models with robust system dynamics during the search process. in the sos case, the parameters indicated in table can be used to derive the scaffold of the sos network, as shown in figure the most sensitive parameters identified by the proposed approach. the letter g indicates that gi,j is the most influential parameter for gene i; while the letter h indicates that hi,j is the most influential parameter. the asterisk means the regulation known in the literature. to compare the gene regulation events identified by our approach to those found in the literature, we highlighted the gene interactions in table and summarized the results in table the results show that, of the thirteen gene regulation events found in our experiments, ten matched regulatory relationships known from other studies. for example, in figure the negative regulation of lexa, uvra, and uvrd by lexa has been successfully identified as the most crucial interaction in determining system dynamics. still, it is notable that the known negative regulation of lexa by reca was not recognized. the main reason for this is that our approach chose sensitive parameters from the kinetic orders of the s-system model on the gene level, but the above relationship was in fact interpreted as a regulation of protein lexa by protein reca. the concentration of protein reca depends on the sensing of dna damage, and the interaction between reca and lexa therefore depends on the events occurring between dna damage and dna repair. because the on or off state of dna damage sensing was not directly described in the s-system model, the regulation of lexa by reca on the protein level was not observable on the basis of kinetic orders. in other words, the absence of this regulation was due to the model representation, rather than the inference approach. this shows that, to infer the sos repair system accurately on both the protein and gene levels, a more comprehensive mathematical model is needed to describe the regulatory details. with only one exception that was caused by limitations of the model itself, our approach has shown its strength in inferring real gene system that show the expected network behavior and internal network robustness. summary of critical gene regulatory relationships obtained by the proposed framework with respect to the results in the literature. the parameters gi,j and hi,j are kinetic orders of the s-system. after identifying the most sensitive parameters of the sos repair system, we again analyzed the gene regulation relationships obtained from the runs to observe the correlation of parameters. the importance of investigating the relationships among parameter values of genes has been emphasized in the study of systems biology, especially when the boundary of a parameter value is taken into consideration. as presented and discussed in the second experimental phase, the qualitative relationships among genes can be regarded as special patterns that determined the system dynamics of a candidate model. by exploiting these patterns, the inference algorithm can refine the values of kinetic orders to obtain better solutions. in the case of the sos repair system, the patterns of regulatory interactions can be derived as in other cases, by recording and analyzing the values of the most sensitive parameters. taking gene lexa as an example, the algorithm categorized genes lexa and uvra as the critical gene regulators for lexa and two regulatory relationships, lexa -| lexa and uvra -| lexa can be established on the basis of this information. figure illustrates the parameter values related to the patterns derived from the twenty runs. in this figure, the x-axis indicates the run number, and the y-axis represents the parameter values of the specific pattern obtained in each run. in this figure, the upper and lower parts present the results of inferring the network without and with using the sa method, respectively. the regulation patterns for gene lexa are summarized in table as shown in table for experimental runs using setting fifteen out of twenty runs consistently produced the same pattern uvra -| lexa > lexa -| lexa. this pattern can guide the inference algorithm to find better solutions. in contrast, the qualitative relationships between the two relationships uvra -| lexa and lexa -| lexa obtained from the runs using setting are inconsistent and unstable. the above results show that the proposed approach can derive parameter bounds to restrict the variation of lexa, and it can therefore infer models with better fitness and more stable system dynamics. without using the sa method, the inference algorithm alone cannot guarantee the preservation of useful patterns during the iterative inference process. pattern-ns1: pattern-ns2: pattern-s1: the patterns of gene regulations found for target gene lexa. tables list the regulation patterns analyzed for four other target genes in the sos system: lexa, uvra, uvrd, and umud. figures show the parameter values involved. the genes reca and polb were ignored because there was only one gene regulatory relationship found for each of them. this means that no qualitative relationship can be established for these genes. the results in tables are consistent with table which again shows the success of using sa in network inference. pattern-ns1: pattern-s1: the patterns of gene regulations found for target gene uvra. pattern-ns1: pattern-s1: pattern-s2: the patterns of gene regulations found for target gene uvrd. pattern-ns1: pattern-ns2: pattern-s1: pattern-s2: the patterns of gene regulations found for target gene umud.  in this work, we developed a sensitivity-based incremental evolution approach to cope with one important issue, network robustness, which has not been addressed in gene network inference. our approach included two parts. the first part was a sensitivity analysis method that was used to select sensitive parameters for deriving value bounds of these parameters. the second part was an evolutionary algorithm that took the bounds as constraints to perform parameter optimization. this process leads to inferred networks that are robust and produce the desired behaviors. to validate the proposed approach, a series of experiments were conducted to evaluate the external behaviors and internal robustness of the networks inferred by different methods. the results show that the proposed sa-based approach outperformed other methods in all datasets. moreover, we analyzed in detail the results obtained from real time-series expression data for the sos repair system. the analyses indicate that our approach can identify the critical parameters and use them to establish regulatory relationships among genes. by enforcing these relationships in the repetitive search process, our approach can successfully infer robust networks. competing interests the authors declare that they have no competing interests. authors' contributions yh undertook the experimental implementation, made statistical analysis and wrote a part of the manuscript. wl conceived the project, designed the algorithm, wrote a part of the manuscript, and made modifications as well as final revisions. 
we have conducted a series of experiments to evaluate the external behaviors and internal robustness of the networks inferred by the proposed approach. our approach includes a sensitivity analysis procedure to iteratively select the most influential network parameters, and it uses a swarm intelligence procedure to perform parameter optimization. we propose an incremental evolution approach for inferring grns that takes network robustness into consideration and can deal with a large number of network parameters. inferring robust networks with desired behaviours remains challenging, however. sensitivity analysis is crucial to identifying the most sensitive parameters that govern the network dynamics. this problem is related to network dynamics but has yet to be investigated using network modeling. the results and analyses have verified the effectiveness of our approach. the experimental results show that the proposed approach can successfully infer robust grns with desired system behaviors. reconstructing gene regulatory networks from expression data is one of the most important challenges in systems biology research. it can further be used to derive constraints for network parameters in the network reconstruction process.
 reconstructing gene regulatory networks from expression data is one of the most important challenges in systems biology research. many computational models and methods have been proposed to automate the process of network reconstruction. inferring robust networks with desired behaviours remains challenging, however. this problem is related to network dynamics but has yet to be investigated using network modeling.  we propose an incremental evolution approach for inferring grns that takes network robustness into consideration and can deal with a large number of network parameters. our approach includes a sensitivity analysis procedure to iteratively select the most influential network parameters, and it uses a swarm intelligence procedure to perform parameter optimization. we have conducted a series of experiments to evaluate the external behaviors and internal robustness of the networks inferred by the proposed approach. the results and analyses have verified the effectiveness of our approach.  sensitivity analysis is crucial to identifying the most sensitive parameters that govern the network dynamics. it can further be used to derive constraints for network parameters in the network reconstruction process. the experimental results show that the proposed approach can successfully infer robust grns with desired system behaviors. 11- august the international conference on intelligent computing zhengzhou, china 