Sensitivity analysis approaches applied to systems biology models

: With the rising application of systems biology, sensitivity analysis methods have been widely applied to study the biological systems, including metabolic networks, signalling pathways and genetic circuits. Sensitivity analysis can provide valuable insights about how robust the biological responses are with respect to the changes of biological parameters and which model inputs are the key factors that affect the model outputs. In addition, sensitivity analysis is valuable for guiding experimental analysis, model reduction and parameter estimation. Local and global sensitivity analysis approaches are the two types of sensitivity analysis that are commonly applied in systems biology. Local sensitivity analysis is a classic method that studies the impact of small perturbations on the model outputs. On the other hand, global sensitivity analysis approaches have been applied to understand how the model outputs are affected by large variations of the model input parameters. In this review, the author introduces the basic concepts of sensitivity analysis approaches applied to systems biology models. Moreover, the author discusses the advantages and disadvantages of different sensitivity analysis methods, how to choose a proper sensitivity analysis approach, the available sensitivity analysis tools for systems biology models and the caveats in the interpretation of sensitivity analysis results.


Introduction
Cell behaviours are not only determined by the characteristics of individual biological components but also by the interactions of such components acting together as a system.Systems biology approach with mathematical models has emerged as a powerful tool in studying the dynamics of the biological systems, because it provides a way to predict emergent network properties and to uncover the principles of cellular networks by merging prior knowledge with experimental data and model simulations [1].The development of predictive dynamic models requires the information about the initial conditions and kinetic parameters that characterise the biological systems.Unfortunately, the parameters of systems biology models are difficult and sometimes even impossible to measure with biological experiments.In addition, some parameter values have large variations among different experimental conditions.Thus, our confidence on the model predictions is limited due to the uncertainties of the model parameters.Sensitivity analysis is a classic technique to determine how the fluctuations in mathematical model outputs can be apportioned to the variations in the model inputs [2].
Sensitivity analysis plays an important role in dynamic analysis of systems biology models [3].First, sensitivity analysis can provide valuable insights about how robust the model outputs are with respect to the changes of model inputs and which model inputs are the key factors that affect the model output.In addition, sensitivity analysis is important in model development.It helps us to test whether a model prediction is dependent on the model assumptions and guide parameter estimation and experimental design [4].Sensitivity analysis is also useful for model simplification by quantifying the dependence of model outputs on its inputs.The values of the insensitive parameters might be fixed and some processes may be simplified or eliminated [5].Last, but not the least, sensitivity analysis can guide experimental analysis.With sensitivity analysis, one can pinpoint which model inputs contribute most to the variation in model outputs (experimental observations).The most sensitive model input parameters and their corresponding biological processes are the potential targets for further experimental analysis.
Sensitivity analysis has a long history and it has been widely applied in different fields such as environmental modelling study [6], economic modelling for decisions making [7], chemical kinetics [8][9][10] and biological modelling analysis [11,12].In this review, we particularly focus on the sensitivity analysis approaches applied to systems biology models.In addition, we limit the discussion of relevant model inputs to the kinetic parameters and the initial conditions of systems biology models.Although the systems biology models discussed here are in the format of ordinary differential equations (ODEs), the sensitivity analysis approaches are generally applicable to other formats of biological models.
There are two types of sensitivity analysis approaches: local and global sensitivity analysis.Local sensitivity analysis is a common approach that the sensitivity of a model output is performed by computing the first-order partial derivatives of the system output with respect to the input parameters, which can be viewed as the gradients around the multidimensional reference parameter space [9].The second type of methods is global sensitivity analysis, which is used to quantify the overall effects of the model inputs on the model output by perturbing model input parameters within large ranges.The introduction of global sensitivity analysis approaches will be presented in the next section.
The systems biology models discussed here is a system of ODE that is dependent on a certain parameter set p and initial conditions y i (0), which is Local sensitivity analysis studies the changes in the model outputs (e.g.ODE variable y i ) with respect to model input (e.g.parameter p) variations around a local point in the parameter space, which are quantified by the sensitivity coefficients.Mathematically, the sensitivity coefficients are the first-order derivatives of model outputs with respect to the model parameters where y i is the ith model output and p is the model input parameter.

Finite difference approximation
There are different methods for the calculation of the derivative.A simple method is computing the derivative by finite difference approximation In reality, the accuracy of numerical approximation depends on perturbation size (Dp) and the precision of the numerical simulation for the model system.In most cases, the ODE models do not have analytic solutions and they are solved with numerical approximation.Therefore the parameter perturbation should be small enough to have a small error in the finite difference approximation, and large enough to reduce the dependence on the simulation inaccuracies from ODE solver.As a rule of thumb in practice, we can start from a small perturbation of the model input (e.g.Dp ¼ 0.001 × p), then change the perturbation size (Dp) smaller or bigger and compare the results.If the results of sensitivity do not change significantly, the local sensitivity analysis result is robust and trustable.

Direct differential method
Another way to compute the derivatives is the direct differential method, which solves the differential equations for the sensitivity coefficients [9,13].These differential equations are differential output y i with respect to the model input parameter p over time.
We can further derive by using the chain rule of differentiation.
It is worth noting that (∂f i /∂y j ) in ( 5) is an element of the Jacobian matrix ( J ) of the ODE model.Thus, (5) can be rewritten as The matrix f p , J and S have the following definition The initial conditions of the ODE system (6) can be determined by the relationship between the model input parameter p and the initial conditions of y variables [9].If model input p is not one of the initial conditions, thus S i (0) ¼ 0 for all i based on the definition with (2).If the model input p is the initial condition for y k , then S i (0) ¼ 0 for i = k and S i (0) ¼ 1 for i ¼ k.The sensitivity of ODE system (6) can be numerically solved with the CVODES solver in the SUNDIALS library [14,15].The advantage of the direct differential method over finite difference approximation approach is that it avoids the difficulty of selecting the variation size Dp by trial and error [9].Furthermore, the sensitivity of different variables with respect to a certain parameter p can be solved simultaneously by the direct differential method.The disadvantage of the direct differential method is that the Jacobian matrix ( J ) needs to be defined and it is timeconsuming especially for large-scale problems.

Adjoint sensitivity analysis
The direct differential method described in the previous section becomes computation-demanding if the sensitivities for many parameters are analysed [15,16].To reduce the computation cost, an alternative adjoint method was developed.In the adjoint method, a Green's or kernel function matrix K(t, t) is constructed for (6), which satisfies dK(t, t) dt = J (t)K(t, t), t ≤ t (8) The sensitivity of S with respect to the parameter p can be rewritten with the constructed Green's function (see [17] for details) In the integral part of (10), the first argument (t) of K is fixed, whereas the second argument (t) is varying.To evaluate the integral part in (10) more efficiently, an adjoint strategy is employed.The differential equations for the adjoint Green's function K * (t, t) are defined as and solved backwards in time from t to 0 using the following identity The numerical calculation of (10) was described in [16,17].The sensitivity coefficients obtained by (10) are the same as those obtained with (6).The difference between the adjoint method and direct method is the numerical methods.The adjoint method is computationally less expensive than direct method when the number of parameters is greater than the number of species (for detailed explanation, see [9]).

Metabolic control analysis (MCA)
Depending on the selected model output and model input, several types of sensitivity coefficients are defined in different sensitivity analysis studies.MCA is a standard mathematical framework for quantifying how the properties of a biochemical reaction network depend on the network parameters.MCA was originally developed to study the metabolic networks [18 -22], but the concept was later applied to other biological networks such as cell signalling, genetic networks and other biological processes [23 -25].
The central idea of MCA analysis is the control coefficient that quantifies the relationship between model output y (e.g.concentrations and fluxes) and input x, which is a normalised derivative and defined as [18 -20, 22, 26] Different control coefficients are defined in MCA analysis based on the choice of model output y and input x (Table 1).
p-elasticity coefficient (parameter-elasticity coefficient) the reaction rate v i the parameter p j [27,28] response coefficient the steady state of state variables or fluxes the parameter p j [26] It is worth noting that different names of control coefficients have been reported during the historical development of MCA theory.Some researchers made agreement on the uniform nomenclatures for principle types of coefficients [32].Two types of coefficients were defined: global and local coefficients.The global coefficients refer to the sensitivity measurements in the entire system.Global coefficients are usually calculated from a given steady state.In contrast, the local coefficients stand for the individual functional entity that is 'isolated' from the system [32,33].The meanings of 'local' and 'global' coefficients defined here are with respect to the system, rather than the parameter space.For instance, elasticity coefficients are local coefficients, while control coefficients and response coefficients are global coefficients [26].

Global sensitivity analysis
The definitions of sensitivity coefficients described above indicate that local sensitivity analysis is estimating the derivatives at a particular point in the model parameter space.However, it is likely that biological model inputs such as rate constants and initial concentrations are varied within a large range in different cell types and cellular environments.For this reason, global sensitivity analysis approaches were applied to quantify the sensitivities of the model outputs with respect to large variations of the model input parameters.This section aims at introducing the principles of different global sensitivity analysis approaches instead of explaining the detailed algorithms for these methods.The examples of global sensitivity analysis applied to systems biology models are described in Table 2.

Parameter space sampling: Latin hypercube sampling
Global sensitivity analysis approaches require taking samples to mimic the variations in model inputs.Therefore parameter space sampling is essential for global sensitivity analysis.Simple random sampling is an easy method to generate samples from a predefined distribution.However, it is shown that simple random sampling may require a large number of samples to cover the entire range of the model inputs [49], which results in the computation-demanding issue for complex models.Latin hypercube sampling is an alternative method because it is an efficient method to sample random parameter vectors while guaranteeing that individual parameter ranges are evenly covered.This method was first developed by McKay Conover in 1979 [50] and later extended by Iman et al. [51].In general, Latin hypercube sampling can generate the input distribution with fewer sampling iterations than the simple Monte-Carlo sampling to achieve a similar accuracy.
The procedure of Latin hypercube sampling for selecting M different values from each of N variables x 1 , x 2 , . .., x n can be summarised as: (i) Divide the range of each variable into M equally intervals; (ii) From each interval, randomly select a value with respect to the probability density in the interval; (iii) The M values thus obtained for x 1 are paired randomly with the M values of x 2 .These M pairs are combined in a random manner with the M values of x 3 to form M triplets, and so on, until a (M × N ) matrix is formed.The implementation of Latin hypercube sampling is available as functions in different programming languages.For example, Latin hypercube sampling can be implemented by using 'lhsdesign' function in MATLAB [52].
The advantage of Latin hypercube sampling compared to simple random sampling is that it can guarantee that the parameters are uniformly sampled from the corresponding distribution.For example, if we generate 10 samples for the two variables, x and y, 10 (x, y) data points are to be sampled.The two variables are assumed uniformly distributed between 0 and 1.If the samples are generated with Latin hypercube sampling method, the sampling points  [43] thermodynamic gene expression model [45] RS-HDMR high genetic circuit model [46] PI3 K/Akt signalling pathway model [47] tumour growth and dendritic cell therapy [48] IET Syst.Biol., 2011, Vol. 5, Iss. 6, pp.336-346 339 doi: 10.1049/iet-syb.2011.0015 of x and y will be equally distributed between 0 and 1 (Fig. 1a), which means that x and y will have one and only one sample at each of the 10 intervals (0, 0.1], (0.1, 0.2], (0.2, 0.3], . .., (0.8, 0.9], (0.9, 1).In contrast, if the samples are drawn from a uniform distribution by simple random sampling, more than one sample may fall into some intervals while no samples might hit other intervals (Fig. 1b).

Multi-parametric sensitivity analysis (MPSA)
MPSA method (also called regionalised sensitivity analysis) was first proposed by Hornberger and Spear [6] in the field of hydrology.It is a sensitivity approach based on Monte-Carlo filtering.Cho et al. and Zi et al. [34,35] first applied this method to study cell signalling pathways.The MPSA method was later used in the analysis of other systems biology models [11,36,37].
The basic principle of MPSA method is to evaluate the model output for the randomly generated parameter sets and then evaluate the parameter sensitivity based on Kolmogorov-Smirnov statistics by classifying the parameter sets.The modeller selects the parameters and randomly generates the parameter space from a given distribution.Latin hypercube sampling method can be used for parameter sets sampling.In the next step, an objective function (model output) is computed for each random parameter set.The random parameter set is then classified as acceptable or unacceptable by comparing its objective function value to a threshold, which can be defined as the average of all the objective function values.If the objective function value is smaller than the threshold, the parameter set is classified as 'acceptable' or 'behavioural'; otherwise it is 'unacceptable' or 'non-behavioural'.For each parameter, the parameter sets are sorted with the increment of this parameter values.Correspondingly the cumulative frequency is calculated for both acceptable and unacceptable cases.Finally, parameter sensitivities are measured by the maximum vertical distance of the two cumulative frequency curves based the Kolmogorov -Smirnov statistics [35].The calculated MPSA sensitivities are between 0 and 1. Larger parameter sensitivity indicates that the corresponding parameter variation has large impact on the defined model output.The detailed algorithm and implementation of the MPSA method have been described before [34,35].

Partial rank correlation coefficient (PRCC) analysis
PRCC analysis is a sensitivity analysis method that calculates the partial rank correlation coefficient for the model inputs (sampled by Latin hypercube sampling method) and outputs [11,38,53,54].The PRCC method assumes a monotonic relationship between the model input parameters and the model outputs.Therefore the PRCC method can be used only when the model outputs are monotonically related to the model input parameters.The monotonic relationship between the input and output can be examined by inputoutput scatter plots [38].
The calculated PRCC values are between 21 and 1 and they are comparable among different model inputs [11].The sign of the PRCC values shows the qualitative relationship between the model input and model output.A positive PRCC value implies that when the corresponding model input increases, the model output will also increase.A negative PRCC value indicates a negative correlation between the model input and output.The magnitude of the PRCC sensitivity measures the importance of the model input in contributing to the model output [38].Details about the PRCC method were described in previous publications [11,38,54].

Morris sensitivity analysis method
The Morris method is a screening sensitivity analysis that is based on an elementary effect [55].The elementary effect is calculated by changing one parameter at a time.The parameter is sampled p discrete times from a predefined distribution.The sensitivity analysis of a model output for k parameters requires a k × p grid.Considering the sensitivity analysis of a model output y for the k model parameters, the ith model parameter can be scaled in the interval [0 1] and may take values from {0, 1/( p 2 1), 2/( p 2 1), . . ., 1}.The elementary effect of the ith parameter x i is defined as where x is the reference parameter set u ¼ [x 1 , x 2 , x k21 , x k ], x i ≤ 1 2 D and D is a predefined multiple of 1/(p − 1).The average of the elementary effects quantifies the importance of the parameters for the model output, while the standard deviation of the elementary effects indicates the non-linear effect of the model parameters on the model output [55,56].A large mean of elementary effects implies that the corresponding model parameter has an important influence on the model output.A high standard deviation of elementary effects indicates that either the parameter is correlated with other parameters or the parameter has nonlinear effects on the output [56].The Morris method has a low computational cost, which is appropriate to study largescale models.

Weighted average of local sensitivities (WALS)
Bentele et al. [4] proposed a WALS coefficient to measure the global parameter sensitivities.The concept of this method is similar to the Morris sensitivity analysis method.In this approach, local sensitivity coefficients are calculated at different random locations within the parameter space.For each parameter, the occurrences of local sensitivity coefficients in different parameter sets are weighted based on a Boltzmann-distribution weighting function, exp(2E/ k b T ), where the variable E can be defined as the weighted least squares error (WLSE) between the perturbed model output and reference model output (or experimental observation).The denominator 'k b T ' is a customisable scaling factor, which can be defined as the minimum of WLSE for all the sampling parameter sets [11,12].The weight factor for the local sensitivity of parameter p in the parameter set x k is defined as The global sensitivity coefficient of parameter p is calculated as where S k p is the local sensitivity of parameter p in the parameter set x k , w k p is the weight factor defined by ( 16), and N is the number of sampling parameter sets.

Sobol sensitivity analysis method
The Sobol method is a variance-based sensitivity analysis approach that makes no assumptions on the relationship between the model inputs and outputs [57].The basis of the Sobol method is the decomposition of the variance of the model output function f (x) into summands of variances in combinations of input parameters in increasing dimensionality [57] The total variance D is defined as The partial variances are computed from each item in (18) and defined as The variance of the model outputs D and D i 1 i 2 ...i s can be approximated by Monte-Carlo integrations.The Sobol global sensitivity indices are calculated by The term S i 1 i 2 ...i s gives the fraction of the total variance which is apportioned to the individual model parameters or the combination of them.For example, S i ¼ D i /D is called the first-order sensitivity index, which quantifies the contribution of the parameter x i to the output variance.Homma and Saltelli suggested the total effect sensitivity index as an extension of the Sobol sensitivity indices, which is defined as the ratio of the sum of the related partial variances [58] S T i = S i + S ci (22) where S ci is the complementary set of practical sensitivity indices that parameter x i appears.For example, S T 1 = S 1 + S 12 + S 13 + S 123 is the total effect index of model parameter x 1 for three-model input parameters.The total effect sensitivity index quantifies the overall effects of a parameter, in combination with any other parameter(s), on the model output.The Sobol method is computation-demanding because it requires a large number of model simulations with the winding stair algorithm [11,12].

Fourier amplitude sensitivity test (FAST)
Another variance-based global sensitivity analysis is FAST.The FAST method was proposed by Cukier and Shuler et al. [59 -61] to study the chemical reaction systems in 1970s.The FAST method assumes that each model parameter is statistically independent from all the others [9].The parameter is sampled from the following transformation function where p 0 i is the reference value for parameter i. w i is an element in a set of linearly independent integer frequencies.s is a scalar variable.G i is a defined transformation function that transforms the probability density of the parameter into s space [62].
According to (24) and the application of Parseval's theorem, the output variance D can be approximated by performing a Fourier analysis as [63] where A j and B j are the Fourier coefficients The partial variance in FAST method is approximated by The FAST sensitivity index is calculated with the following approximation The implementations of the FAST method are complicated and they are available in [54,59,60,62].The computational cost for the FAST method is very high because it needs a large number of model evaluations.
An extended FAST (eFAST) method was proposed by Saltelli et al. [64], in which a transformation function G i was recommended as The eFAST method has a better transformation function than the classic FAST method because it provides uniformly distributed samples for the parameters [64].

Random sampling high-dimensional model representation (RS-HDMR)
HDMR is another variance-based method, which was introduced for improving the analysis of the input -output behaviour in high-dimensional systems [65,66].The HDMR method has been applied to the global sensitivity and uncertainty analysis of the mathematical models.It expresses the model output f (x) as a combination of the model inputs x, which is a statistical ANOVA decomposition [67,68] where the component function f 0 is a constant corresponding to the zero order effect and the first-order component function f i (x i ) describes the independent contribution of the ith parameter x i to model output f (x), the second-order component function f ij (x i , x j ) represents the pair interactive contribution of the two model parameters x i and x j to model output f (x) and the last term gives the cooperative contributions of all the input variables to the model output.
Previous work indicates that HDMR expansions to the second order can provide a satisfactory description of the model output for many systems [65,67].
The sensitivity analysis method by RS-HDMR is based on the RS-HDMR component functions, which have the following forms where x i , x ij are x without the element x i and x i , x j .dx i and dx ij represent the product dx 1 dx 2 . ..dx n without dx i and dx i dx j , respectively.For RS-HDMR component functions, the model input x i is rescaled such that 0 ≤ x i ≤ 1.
Using the orthogonality of the component functions, RS-HDMR defines the total variance s 2 [47, 65] The RS-HDMR sensitivity indices, S i (i ¼ 1, 2, . . ., n p ), are defined as the portion of the total variance s 2 represented by the variance of the ith component function.In addition, total sensitivity indices S T i can be defined as the sum of the sensitivity indices S i describing the first and the secondorder (or higher orders if wanted) component functions of a model input x i .The total sensitivity indices S T i describe both independent and higher-order effects of the model input x i on the model output [47].The approaches for calculating RS-HDMR component function were described in the previous publications [67, 69 -71].
The advantage of RS-HDMR global analysis method is that it works well even when model input parameters contain large uncertainties, which is critical for the application to systems biology models where model parameters are not precisely identifiable [46].

Sensitivity and identifiability of biological models
In this review, we have introduced different sensitivity analysis methods that are applied to systems biology models.Biological models seem to be very robust to the changes of many parameters.Gutenkunst et al. investigated the sensitivities of 17 published systems biology models and studied the model output variations to the parameter changes.They found that systems biology models exhibit sloppy sensitivity spectra because the behaviour of the models is very insensitive to many parameter changes [72].The general sloppiness in biological model parameters implies that it is difficult to uniquely determine the model parameters by fitting to a few experimental data.However, the parameters of some biological models can be calibrated and well-constrained with multiple experimental data sets.For these data-based models, local sensitivity analysis may be very powerful to pinpoint the key factors that affect the model behaviours.The sensitive parameters might be useful for further investigation of the biological mechanisms.For example, local sensitivity analysis has recently been employed to identify the critical signalling steps that are responsible for the linear detection of the Epo ligand concentrations [73] and the ultrasensitivity of signalling responses in TGF-b pathway [74].Therefore the precision of model prediction from local sensitivity analysis is dependent on the identifiability of model parameters.On the other hand, local sensitivity analysis can be used to construct Fisher information matrix (FIM) for parameter identifiability analysis, which determines whether certain parameters can be identified from the experimental data [5,75,76].The FIM is a function of the local sensitivity matrix.The 95% confidence intervals of the estimated parameters can be approximately calculated from the Crame ´r -Rao lower bound for the variance of the estimators based on the FIM (see [76] for details).

Pros and cons of global sensitivity analysis approaches
Here, we discuss the advantages and limitations of some global sensitivity methods as a guidance for the application of sensitivity analysis in biological models.

MPSA method:
The MPSA method is easy to implement.It has some global properties: the whole range of input parameter values is covered and all the input parameters are varied at the same time.The Kolmogorov -Smirnov analysis in the MPSA method relates not only to the main effects as variance-based methods, but also highlights certain types of interaction effects [2].While the MPSA method can be very useful in providing the overall impact of the input parameters on the model output, it also has some limitations: (i) It can be subjective as the modeller need to choose a threshold for defining the 'acceptable behaviour'.(ii) It is not always obvious to select a proper distribution of parameters to guarantee the Monte-Carlo simulations generate sufficient 'acceptable' hits that make the statistical analysis meaningful [77].Previous practice with hydrological models has shown that the fraction of 'acceptable' behaviours is barely larger than 5% over the total simulations for large models (e.g. with number of input parameters .20),indicating a lack of statistical power [2,78].(iii) If some parameters are highly correlated with each other, the marginal distributions cannot separate under the acceptable classification and thus underestimate the importance of some parameters.The Kolmogorov -Smirnov test is sufficient to ensure whether a parameter under analysis is important.However, it does not provide a necessary condition for importance.In other words, the unimportant parameters identified by the MPSA method do not make sure that they are trivial to the model output.Such feature of the MPSA method implies that it may perform incomplete assessment of the analysed parameters and might lead to some false-negative results.Further complementary analysis (e.g. by applying other global sensitivity analysis approaches) is necessary to verify that the unimportant parameters are not involved in the interactions [2].

Morris method:
The main advantage of the Morris method is its computational efficiency due to the relatively low computational cost.A drawback of the Morris method is that the sensitivity measure is only qualitative as it only gives an overall measure of the interactions.When the model is non-monotonic, the distribution of the elementary effects, d(x i ), may have negative elements, some effects may cancel each other out when computing the mean [79].Thus, the mean measure is not reliable for ranking the importance of the parameters.In the implementation of the Morris method, the choice of the size of the levels p and the parameter D may also affect the robustness of the results.

PRCC method:
The PRCC method is more robust than simple correlation coefficients analysis approaches as it uses rank transformation statistic.Rank transformation of the data can transform a non-linear but monotonic relationship to a linear relationship.This increases the applicability of PRCC sensitivity analysis method as linear relationship is not required for it.However, the performance of the PRCC method is limited by another assumption that the model output should be monotonically related to the model input parameters.Thus, the PRCC analysis results may be misleading or not accurate in case that nonmonotonicities exist in systems biology models [54].

4.2.4
Variance-based sensitivity methods: Saltelli et al. [2] suggested to use the variance-based sensitivity analysis approaches (e.g. the Sobol, FAST, HDMR methods) when such application is possible.This recommendation is based on the advantages of the variance-based methods: independence on model linearity or monotonicity, capability to obtain the impact of the full range of each input parameter variation and allowing the interaction effects among input parameters.The main drawback of the variance-based approaches is their high computational cost because they require more model evaluations than other types of methods.It can become prohibitive for the application of variancebased methods in the case of computationally laborious models including large amount of parameters and state variables.

Choice of sensitivity analysis methods
A general question for applying sensitivity analysis to systems biology models is: how to choose a sensitivity analysis method?The answer is dependent on many factors, for example, the number of input parameters, the computational cost of the running model, the uncertainty of the model parameters, the correlation and interactions among the input parameters and the features of the model (e.g.linearity or monotonicity).There is no absolute good sensitivity analysis method for all the models because most of the methods have their pros and cons as mentioned in the previous section.Here, we list some guide rules for the choice of sensitivity analysis approaches (for more details, see [2,80]).

Local or global sensitivity analysis?:
In general, local sensitivity analysis is relevant to study how small perturbations of the model inputs affect the model behaviours, while global sensitivity analysis is more favourable to investigate the impact of large model input variations.The derivative-based local sensitivity approaches are attractive due to their low cost in computation time.Yet, the fatal limitation of a local sensitivity approach is that it is unwarranted when the model parameters are nonidentifiable or uncertain (e.g. with uncertainties of different order of magnitude).In other words, the derivatives are only informative at the base point in the reference parameter set and do not explore other positions in the parameter space [2].Local sensitivity analysis is proper for linear systems or the analysis of the input parameters that are identifiable.As the predictability of global sensitivity analysis does not rely on the identifiability of model parameter set, global sensitivity analysis approaches are more suitable for determining which of the uncertain input parameters are more important for the variation in the output of interest.

Choice of global sensitivity analysis methods:
First, the computational cost of the sensitivity analysis method should be taken into account.For example, the variance-based methods require high computational cost.Therefore these methods might not be a good choice for large models with many parameters to be analysed.The screening-based Morris method is useful as a first step in dealing with computation-demanding models containing a large number of input parameters.In addition, the features of the model should also be considered.As a starting point, one can judge the non-linearities, non-monotonicity and correlation between the input -output parameters by coefficient of determination (R 2 -value) or scatter plot.For linear models, all the global sensitivity analysis methods are applicable.Sampling-based approaches (e.g.MPSA, WALS and PRCC) are recommended because they have low computational cost.For non-linear but monotonic models, the PRCC method can perform well.For non-monotonic models, the Morris method will produce poor results because the importance of the parameters might be underestimated due to the cancelling out in computing the mean of the elementary effects.In addition, the PRCC method should be avoided because it assumes a monotonic model input -output relationship.In this case, the variancebased approaches such as the Sobol method, FAST and HDMR are recommended if the application of these methods is possible [80].

Application of Latin hypercube sampling method:
Latin hypercube sampling is useful for the sampling-based sensitivity analysis methods as it ensures the samples to be evenly drawn from the full range of the desired distribution functions.A small number of samples are needed to mimic the distribution functions and correspondingly the number of required model evaluations can be reduced.A drawback of Latin hypercube sampling is that it can yield a biased estimation of the variance of the distributions.Thus, Latin hypercube sampling should be avoided for generating sample distributions in the variancebased sensitivity analysis methods [80].

Timing matters for sensitivity analysis
As the impact of model parameters on the model output changes over time, time-dependent parameter sensitivity analysis has been proposed to study the effect of parameter variation on model output at different time [29][30][31].A parameter may have positive impact on the change of model output at early stage, but its effect can switch from positive to negative due to the complex feedbacks in the biological network.Therefore one needs to know not only which parameters are critical for affecting a model output, but also at which time point they are important.Moreover, traditional sensitivity analysis methods are computed for persistent parameter change, which means that parameters are constantly changed over time.Based on the consideration of time effect, Perumal and Gunawan recently pointed out that the sensitivity of a parameter is also dependent on the time when the parameter is changed.Accordingly, Perumal and Gunawan proposed a new method called impulse parametric sensitivity analysis, which is useful to understand when the parameter variation will cause the change of systems dynamics [81].

Software tools for sensitivity analysis of systems biology models
The numerical implementations of sensitivity analysis methods are tedious especially for some global sensitivity analysis algorithm.This imposes a demand for systems biology community to develop tools that allow automatic sensitivity analysis of systems biology models.Several software tools have been developed to perform sensitivity analysis for biological models, for example, BioSens [82], COPASI [83], SBML-SAT [12], Systems Biology Toolbox 2 [84], SensSB [85], TinkerCell [86] and Virtual Cell [87].These software tools can perform local sensitivity analysis and some of them are able to execute global sensitivity analysis (e.g.SBML-SAT, SensSB and Systems Biology Toolbox 2).Recently, a new method and tool have been developed to map the sensitivity of the model outputs by sensitivity heat maps and a global summation law [88].Most of these sensitivity analysis tools support the models encoded in systems biology markup language (SBML) [89].SBML is a standard markup language for the systems biology community to represent and exchange biological models.The feature of SBML support in the software tools allows the users to perform sensitivity analysis conveniently.Although a few sensitivity analysis tools have been developed, there is no easy-to-use library for the implementation of all the available sensitivity analysis methods.A collection of benchmark examples of systems biology models is still missing for testing and comparing different sensitivity analysis methods and tools.

Final remarks
To sum up, different sensitivity analysis approaches are available for the dynamic analysis of systems biology models.In practice, one needs to carefully choose the type of sensitivity analysis method based on the model features and the purpose of the study.In principle, sensitivity analysis provides a good starting point to identify the input parameters that have strong impact on the model output behaviours.However, sensitivity analysis does not provide a direct explanation or mechanism for such effects [80].One should be cautious about the precision of sensitivity analysis results especially when the model parameters are not identifiable.The sensitivity analysis results should be carefully interpreted.For example, a local parameter sensitivity analysis result might be over-interpreted for predicting the knockout effect of the corresponding biochemical reactions because local sensitivity is based on small parameter perturbations.In practice, perturbation experiments are necessary to test the sensitivity analysis results.Further experimental analyses are important to verify model predictions and understand the dynamics of cellular networks.

Fig. 1
Fig. 1 Sample's design for two variables with different sampling methods a Samples taken with Latin hypercube sampling method.They are uniformed distributed in different intervals for x and y b Samples taken with simple random sampling method.There are no samples falling into the intervals (0.5, 0.9] for variable y IET Syst.Biol., 2011, Vol. 5, Iss. 6, pp.336-346 343 doi: 10.1049/iet-syb.2011.0015& The Institution of Engineering and Technology 2011 www.ietdl.org

Table 1
Definition of control coefficients in MCA analysis

Table 2
Global sensitivity analysis methods applied to