## Abstract

Sensitivity analysis addresses the manner in which model behaviour depends on model parametrization. Global sensitivity analysis makes use of statistical tools to address system behaviour over a wide range of operating conditions, whereas local sensitivity analysis focuses attention on a specific set of nominal parameter values. This narrow focus allows a complete analytical treatment and straightforward interpretation in the local case. Sensitivity analysis is a valuable tool for model construction and interpretation, and can be applied in medicine and biotechnology to predict the effect of interventions.

## Introduction

The models discussed in this *Essays in Biochemistry* volume specify system features by employing a number of parameters, such as kinetic constants, decay rates and association constants. These parameters are treated as constants since they are assumed not to change their values over the timescale of interest. However, these values may be poorly known and in any case are typically dependent on a number of factors which are external to the system of interest and so could vary with experimental, environmental or intracellular conditions. Sensitivity analysis addresses the manner in which model behaviour depends on these parameter values [1–3].

Studies of the effect of parameter variation fall into two categories. Global sensitivity analysis addresses model behaviour over a wide range of parameter values. These ranges are typically based on rough estimates of parameter values, and might be specified by upper and lower bounds, or by probability densities. Statistical methods are often used to guide the sampling of values from within the specified domains of parameter space. An introduction to global sensitivity analysis can be found in [4].

Local sensitivity analysis provides a complement to the broad search procedure adopted in global sensitivity analysis by concentrating attention near a particular point in the parameter space, a nominal operating condition. When addressing model behaviour for parameters near specific nominal values one can safely assume linear dependence on the parameters. This assumption greatly simplifies the analysis and interpretation of the results. There is a long tradition of the use of local sensitivity analysis in biochemistry in work on MCA (Metabolic Control Analysis) [5,6] and BST (Biochemical Systems Theory) [7,8]. The present chapter consists of an introduction to the derivation and interpretation of local sensitivity analysis followed by a brief review of applications of sensitivity analysis in systems biology. For a thorough treatment of global sensitivity analysis the reader is referred to [4]. (See also [9] for a concise introduction.)

## Local sensitivity analysis: derivation

Consider the simple reaction scheme involving a single species *S*:
which might form part of a metabolic chain. Suppose the rate of the first reaction is maintained (by external factors) at a constant level *v*_{1} = *V*, whereas the second follows Michaelis–Menten kinetics:
where *s* = [*S*]. This system is in steady state when *v*_{1} = *v*_{2}, so the steady-state concentration *s*^{ss} satisfies
which can be simplified to
Now, suppose that *V* and *K _{m}* are held fixed while

*V*

_{max}is varied (e.g. through genetic manipulation of enzyme abundance). We then think of this equation as defining the steady-state concentration as a function of

*V*

_{max}, that is

*s*

^{ss}=

*s*

^{ss}(

*V*

_{max}). For the sake of concreteness, take specific values for the other constants, say

*V*= 2 mM/min and

*K*= 1.5 mM. Then the steady-state concentration is given by eqn (1), which is represented by the solid curve in Figure 1.

_{m}This curve, referred to as a continuation diagram, captures the manner in which the steady-state concentration varies as the parameter *V*_{max} takes on different values. Note that the graph does not extend below *V*_{max} = 2 mM/min, as below that value no steady state is possible (since consumption of *S* cannot keep pace with production). This continuation diagram represents a global sensitivity analysis.

Unfortunately, such a complete analysis is rarely applicable to systems of interest. This dependence is generally not directly computable, and so can only be described after simulation of the model at each point along the curve. Moreover this visual representation cannot simultaneously display the system dependence on multiple parameters (beyond two), and it is difficult to compare these dependencies across different parameters.

These issues, which must be faced when attempting any global sensitivity analysis, can be resolved through describing the parameter dependence by a simpler measure, the local parametric sensitivity. As a single number, the local sensitivity coefficient lends itself easily to comparison and representation, and is directly computable as will be shown below. The significant drawback is that it is a measure only of local behaviour, so describes behaviour for parameter values only close to a nominal value.

The local sensitivity coefficient of a steady state *s*^{ss} with respect to a variable *p* is defined as the rate of change of *s*^{ss}(*p*) with respect to *p*, that is, as . This is the slope of the tangent to the continuation curve and indicates the resulting change in *s*^{ss} when a small change Δ*p* is made to *p*, since
In the example above, this sensitivity is given by
In Figure 1, focusing on a nominal value of *V*_{max} = 4 mM/min (which might represent, e.g. the wild-type value), we have a sensitivity coefficient of min, the slope of the dashed tangent line.

In this case, an increase of Δ mM/min in *V*_{max} leads to a 0.75 Δ mM decrease in *s*^{ss}, for *V*_{max} near 4 mM/min. In Figure 1, this corresponds to approximating the solid curve by the dashed line: a valid approximation for values of *V*_{max} near 4 mM/min.

In this example the sensitivity coefficient has units of min, which might seem inappropriate, and is indeed an indication that this is not an ideal measure. In particular, we cannot directly compare this sensitivity to that with respect to *K _{m}*, since is dimensionless. Moreover we can observe that this sensitivity coefficient must be carefully interpreted in the context of the nominal steady state: a change of 0.75 mM will be significant if

*s*

^{ss}= 2 mM, but much less so if

*s*

^{ss}= 200 mM. The root of these issues is that we have been addressing absolute sensitivities, when relative sensitivities, which are computed as ratios of relative changes, are typically more relevant. We scale to yield the dimensionless relative sensitivity coefficient These measurements are sometimes referred to as logarithmic sensitivities or logarithmic gains since In this case we find indicating that a Δ% increase in

*V*

_{max}leads to a 2Δ% decrease in

*s*

^{ss}, for values of

*V*

_{max}near 4 mM/min.

In this simple example we were able to find the sensitivity coefficient after solving explicitly for the steady state in eqn (1). In most models of interest, this explicit solution is not available. In such cases, sensitivity coefficients can be found through implicit differentiation as the following example illustrates.

Consider a simple self-regulating gene circuit, as shown in Figure 2. The gene, which expresses protein *A*, is under the control of a promoter which is inhibited by its own product. A simple model of this system is
where *A*(*t*) is the concentration (nM) of protein *A* at time *t* (min), *d _{A}* (min

^{−1}) describes the rate of degradation of

*A*,

*g*(nM/min) is the maximal rate of expression, whereas

*n*and

*K*(nM) characterize the inhibition of the promoter by

*A*. (This simple model results from abstracting away a number of system features, e.g. mRNA and binding kinetics.)

The steady-state concentration of *A* satisfies
Although we cannot solve this equation for *A* (except in special cases, e.g. *n* = 1), we can find the sensitivity of the steady-state concentration to the model parameters by implicit differentiation. For example, to address the dependence of *A* on *g*, we think of *A* as a function of *g* and differentiate:
Solving for the absolute sensitivity gives
which can be normalized to yield the relative sensitivity coefficient
This formula is rather unwieldy, but allows for easy computation of this local sensitivity once the nominal parameter values are chosen and the steady-state value of *A* is computed. Even without such calculation, one can deduce features of the system directly from the formula. For example, the sensitivity will be positive no matter what (positive) values are chosen for the parameters. Thus the steady-state concentration of *A* increases with increasing *g*, regardless of the operating conditions. This formula can also provide insight into particular operating regimes. For instance, if the steady-state value of *A* is much larger than *K* (i.e. the promoter binding is saturated) so that *A ^{n}* +

*K*

^{n}*A*, then one can verify from eqn (2) that indicating the role of the Hill coefficient

^{n}*n*in setting the sensitivity to

*g*in this case.

For specific values of the parameters we can determine the steady-state concentration of *A* (numerically) and so quantify this local dependence on *g*. For example if *g* = 1 nM/min, *K* = 1.2 nM, *n* = 2 and *d _{A}* = 0.8 min

^{−1}, the steady-state concentration is

*A*= 0.84 nM and . To verify our analysis of the limiting case, we find that if we change

*d*to 0.01 min

_{A}^{−1}, then

*A*= 5.15 nM and so

*A*= 26.5≫1.4 =

^{n}*K*. For these values the sensitivity coefficient is , which is a good fit to our estimate of .

^{n}A more complete model of the activity of this gene includes an explicit description of the mRNA dynamics:
where the mRNA *m* degrades at rate *d _{m}* (min

^{−1}) and is translated at rate

*k*(min

_{A}^{−1}). Again, the steady-state concentrations can be found by setting these rates of change to zero, and sensitivities can be found by implicit differentiation. In this case, and for larger and more complex models, the calculus and algebra required for the sensitivity calculation are cumbersome. Matrix notation can be used to simplify the computation.

The model can be written in vector form as
where **s** is the vector of species concentrations, **f** is a vector function describing the dynamics, and *p* is the parameter of interest. At steady state, we think of **s** as a function of *p*, and arrive at
Differentiating with respect to *p* yields
Solving for the absolute sensitivity yields
provided the matrix , known as the system Jacobian, is invertible (which is the case for most systems of interest). The relative sensitivity coefficients can be found by the appropriate normalization.

In the example above this calculation yields
Fixing parameter values *d _{m}* = 2 min

^{−1},

*k*= 1.5 min

_{A}^{−1}and the other parameters as before, the steady-state concentrations are

*m*= 0.37 nM,

*A*= 0.70 nM and the absolute sensitivities with respect to

*g*are: This analysis can easily be extended to address a number of parameters simultaneously by substituting a vector

**p**of parameters into the formula above. The result of such an analysis is shown in Figure 3, from which one can immediately compare the relative sensitivities of each species to all of the model parameters.

The analytic derivation described above can be readily approximated by numerical means. By simulating the system first at parameter value *p* and then at *p* + Δ*p* one can determine the corresponding steady states and then approximate:
This computational approach is easy to implement and typically gives a good approximation to the sensitivity coefficient, provided Δ*p* is sufficiently small, e.g. 1%. However, one should be aware of the potential for propagation of round-off error in determining ratios of near-zero numbers.

This numerical procedure allows a generalization of the local approach: by choosing larger variations of Δ*p* in eqn (5) one can explore behaviour over a range of parameter values. Such an investigation can complement standard local sensitivity analysis but falls short of a true global analysis since variations are considered one parameter at a time, resulting in a very narrow coverage of parameter space.

Local sensitivities can also be calculated for variables that are dependent on the species concentrations, such as concentration ratios, growth rates, metabolic fluxes or other rates of transmission. Any such variable can be written as *y* = *y*(**s**) and by the chain rule:
This approach will be followed in the next section.

## Local sensitivity analysis: interpretation

To interpret local sensitivity analysis, we take an illustrative example from metabolism. Consider the allosterically regulated metabolic chain shown in Figure 4, where *S* is a boundary species whose concentration is held fixed. Suppose the reaction rates are given as functions of the species concentrations as
with *k*_{01} = *k*_{12} = *k*_{23} = *k*_{34} = 2 nM^{−1}·min^{−1}, *k*_{10} = *k*_{21} = *k*_{32} = 0.1 nM^{−1}·min^{−1}, *K*_{m01} = *K*_{m10} = *K*_{m12} = *K*_{m21} = *K*_{m23} = *K*_{m32} = *K*_{m34} = 1 mM, *S* = 0.1 mM, *e*_{0} = *e*_{1} = *e*_{3} = 6 nM and *e*_{2} = 1 nM. To begin our analysis, consider the case where the allosteric regulation is absent: *p* = 0 mM^{−1}. In addressing metabolic pathways, one is often interested in the steady-state flux through the network which, in this case, is 0.85 mM/min. Using parametric sensitivity analysis we can investigate the dependence of the steady-state flux on the enzyme concentrations. Noting that the other enzymes are six times more abundant than enzyme *e*_{2}, one might presume that reaction *v*_{2} is in some sense the rate-determining step: an increase in *e*_{2} might be expected to have a much more significant effect on the flux than a corresponding increase in the other enzyme concentrations. This is not the case. Indeed the relative dependencies of the steady-state flux *J* on the enzyme abundances are:
(That is, a 1% increase in *e*_{1} produces a 0.69% increase in *J*, and similarly for the other enzymes.) Under these operating conditions the flux is most sensitive to changes in enzyme *e*_{0}, despite the fact that it is far more abundant than *e*_{2}. This non-intuitive result illustrates the power of sensitivity analysis in revealing system behaviour. Even in a simple network such as this, the effect of a perturbation in one reaction reverberates throughout the system until a new steady state is reached. This systemic response can be difficult to predict intuitively in all but the simplest cases.

One must keep in mind that the sensitivity coefficients in eqn (6) are only valid near the given parameter values. For example, if *S* is increased 10-fold to 1 mM, then the sensitivity to *e*_{2} grows to be almost ten times larger than that of the other enzymes, indicating that the rate-limiting intuition is valid under certain operating conditions.

Next, consider the effect of the allosteric inhibition indicated in Figure 4. Setting *p* = 100 mM^{−1} the flux is reduced to 0.29 mM/min. The sensitivities of the flux become:
Without this analysis it would be difficult to work out the effect of introducing this allosteric regulation: a reduction in sensitivity to *e*_{0}, *e*_{1} and *e*_{2} and an increased sensitivity to *e*_{3}. This illustrates a general trend observed by Kacser and Burns [10]: allosteric feedback inhibition has the effect of drawing sensitivity away from those enzymes involved in the ‘loop’ and re-distributing it to those enzymes downstream of the regulating metabolite.

This discussion of ‘re-distributing sensitivity’ can be made precise. The astute reader may have observed that in both of the above cases (eqns 6 and 7) the reported sensitivity coefficients sum to one. That is not a coincidence, but rather follows from the Summation Theorem of MCA, which guarantees this result and so allows direct comparison of sensitivities across different network architectures.

Within MCA, the sensitivities discussed above are referred to as flux control coefficients, a terminology which stresses their role in describing the dependence of pathway flux on individual reaction rates. (To be precise, the sensitivity coefficients derived above are referred to as flux *response* coefficients. The control coefficients can be reached by an appropriate scaling of the response coefficients. In this simple example the scaling factor is one and so the two sensitivity measures coincide.) MCA stresses the relationship (and distinctions) between these system sensitivities and the isolated dependencies of each reaction rate, i.e. the partial derivatives and , which are referred to as elasticities. MCA highlights the manner in which system sensitivities (i.e. control coefficients) can be inferred from elasticities, as was implicit in our earlier derivation of sensitivities in eqn (4). Moreover, MCA stresses the fact that these sensitivity calculations can be carried out in the absence of a complete kinetic description. All that is needed are the elasticities: local descriptions of behaviour which are far more accessible than detailed kinetics.

Thorough treatments of MCA can be found in [5,6]. A complementary treatment of local sensitivities has been applied in the S-system framework of BST [7,8].

## Generalizations of local sensitivity analysis

Local sensitivity analysis as derived in the present chapter describes the effect of small parameter changes on steady-state behaviour. Generalizations of this approach broaden its scope while still addressing small perturbations.

A natural extension is to relax the assumption of steady-state nominal behaviour. Indeed, returning to eqn (3), if we differentiate with respect to the parameter *p* and reverse the order of differentiation, we arrive at
which generalizes eqn (4) and provides a description of the time-varying sensitivity . This approach is discussed in the context of MCA in [11] and in the S-system framework in [12]. The sensitivity coefficients which result from solving the linear differential equation (eqn 8) describe how the transient behaviour of the system is affected when the value of parameter *p* is altered. This analysis is well suited to the study of networks whose dynamics are of prime interest, e.g. signal transduction networks. A recent application to the MAPK (mitogen-activated protein kinase) pathway was presented by Hu and Yuan [13].

As an example, recall the simple genetic circuit in Figure 2. Consider an experiment in which the inhibition has been blocked by a chemical agent and is restored at time zero. Figure 5(a) shows the time-response of mRNA (*m*) and protein (*A*) as their concentrations drop to lower steady-state values under the effect of the inhibition (parameter values as before). Figure 5(b) shows the sensitivity of the two concentrations to a change in parameter *g* over the time course of the experiment.

Analytical treatment of time-varying sensitivities can be demanding, and it is often more efficient to appeal directly to a numerical procedure which compares behaviour at differing parameter values, as in eqn (5). Using parameter variation anywhere from 1% to 200%, this approach has been successfully employed to address sensitivities of the amplitude and period of oscillatory systems (e.g. [14]) and of area-under-the-curve (i.e. integral) measures of both specific variables [15] and of overall system behaviour [16].

As a complement to this non-steady-state analysis, one can address the response of a system to time-varying perturbations. Such an analysis can be carried out through the use of the system frequency response (i.e. via Fourier decomposition). This approach in the context of MCA is described in [17]. An expanded Fourier approach was presented by Liebermeister [18], in which this method is applied to address stochastic parameter fluctuations, i.e. noise.

The discussion in the present chapter has addressed deterministic differential equation-based models. It is becoming increasingly clear that stochastic modelling (as presented in Chapter 3 of this volume) is an essential tool for addressing the behaviour of cellular networks. There is a pressing need for new analytical methods which will aid in the construction and interpretation of these models. The relationship between sensitivities of deterministic and stochastic models is explored by Scott et al. [19]. A direct sensitivity analysis of the chemical master equation underlying most stochastic modelling approaches was developed by Gunawan et al. [20].

The use of local sensitivity analysis is inherently restricted to small regions of parameter space. Although global sensitivity analysis typically makes use of a very different set of analytic tools (primarily statistical), some efforts have been made to marry the two approaches. One such approach is a ‘distributed’ local sensitivity analysis that compares the results of a series of local sensitivity analyses carried out at points sampled from a wide range of parameter space, as in [16,21].

## Sensitivity analysis in systems biology

### Predicting the effect of interventions

Sensitivity analysis plays a valuable role in predicting the effect of perturbations. In fields such as drug design and metabolic engineering, one is interested in achieving a specific change in system behaviour. Local sensitivities provide a first approximation of the effect of various perturbations, and so furnish an initial assessment as to which interventions might be most useful.

For example, Bakker et al. [22] addressed the problem of arresting glycolytic flux in the parasite *Trypanosoma brucei*, which causes African sleeping sickness in humans and nagana in livestock. The metabolism of these protazoa differs sufficiently from their mammalian hosts so that drugs which target their glycolytic enzymes can be used with minimal side effects. Bakker et al. [22] applied a sensitivity analysis (via MCA) to a model of *Trypanosoma* glycolysis and identified enzymes which would make viable targets for pathway inhibition and others whose inhibition might have little effect on the survival of the parasites.

Another example is the work of Comin-Anduix et al. [23] which addresses the dependence of the growth rate of cancerous tumour cells on thiamine. Building on their previous work showing the important role of transketolase in the metabolism of these cells, these authors report an experimentally determined sensitivity coefficient indicating a strong dependence of growth on thiamine availability.

Metabolic engineers are also faced with the question of which perturbation should be made to achieve a specific goal. In this case the systems, typically microbes, are more accessible to investigation and intervention. The goal is often to increase expression of a valuable compound by altering the cellular behaviour through genetic or environmental means. Sensitivity analysis can play the same role here that it does in studies of drug targets. Gárdonyi et al. [24] addressed xylose fermentation by recombinant strains of *Saccharomyces* *cerevisiae*, a necessary process in this organism’ production of fuel ethanol from pentose-rich sources. The authors investigated the sensitivity of the xylose utilization rate to the xylose transport step and were able to experimentally confirm conditions under which two recombinant strains exhibited high or low flux control coefficients for this step.

Hoefnagel et al. [25] consider the production of diacetyl (a flavourant) by the milk-fermenting bacterium *Lactococcus lactis*. The authors used sensitivity analysis together with kinetic modelling to design an engineering strategy which far outperformed previous production methods. Their quantitative analysis revealed that the optimal perturbations were not those suggested by an intuitive approach, and in fact involved perturbing two enzymes which are not part of the main pathway of interest.

The use of MCA in drug development is reviewed in [26], whereas applications to metabolic engineering are discussed in [27]. For a survey of recent uses of MCA in biotechnology and medicine see [28].

### Model construction and analysis

Sensitivity analysis is a standard tool for analysis and design in engineering systems theory [3]. Its application to biochemical models was recognized early on [29]. In recent years, sensitivity analysis has been used for a number of purposes in addressing models of biochemical and genetic networks. Sensitivity analysis can be used to elucidate the behaviour of models. Ihekwaba et al. [14] address a model of the NF-κB (nuclear factor κB) pathway. The authors used sensitivity analysis to identify a small subset of the model parameters which hold the most influence. These parameters characterize the behaviour of just two species in the model, suggesting that a small subsystem may reflect the behaviour of the entire network. A similar analysis was carried out by Stelling et al. [16], in which a distributed local sensitivity analysis is illustrated as a means of investigating the role of network structure in generating system behaviour. A circadian oscillator is used as a test case.

These studies can yield valuable biological insight by revealing the key workings of complex biochemical networks. In addition, they can be used as analytical tools for model reduction, which is often a necessary step in working with large systems. Degenring et al. [30] and Liu et al. [31] used principal component analysis to rank-order parameters by their sensitivities in order to identify model components which could be eliminated with minimal impact on model behaviour.

Sensitivity analysis can also be used as a tool for model calibration: determining numerical values for model parameters. This is a key issue in constructing models with many parameters. Zak et al. [32] made use of sensitivity coefficients to address model identifiability: the question of whether or not parameter values for a model can be determined from a given data set (and more importantly how difficult those determinations would be). Cho et al. [21] take a complementary approach to the same problem. They employ a distributed local sensitivity analysis to arrive at an experimental design which aims to generate data suited for parameter estimation.

In systems biology, models are often used to investigate the robustness of cellular phenomena. Robustness is a rather slippery term which indicates continued performance in the face of an uncertain environment, exemplified by the concept of homoeostasis. Robustness to parameter variation can be equated with low sensitivities. Thus sensitivities can be used for model discrimination, since it is often the more robust model which more faithfully reflects experimental findings, as discussed by Morohashi et al. [33]. A recent example can be found in [34] in which Eiβing et al. use a global sensitivity analysis to compare two models of an apoptotic signalling pathway.

### Continuing developments

Sensitivity analysis is being employed to aid in a novel form of biological design: the construction of synthetic genetic circuits. Feng et al. [35] demonstrated the use of a global sensitivity analysis in identifying components of a circuit design whose optimization would most efficiently lead to improved performance. This sort of analysis, which parallels the use of sensitivity analysis in systems engineering, will no doubt grow in importance as the field of synthetic biology continues to expand.

Another area ripe for future activity is the continued application of global sensitivity analysis to biochemical models. As models continue to grow in complexity, and while high-throughput data continue to be of poor quality for model calibration, these statistical approaches may be our only means for exploring landscapes of high-dimensional parameter spaces. A valuable step in this direction is the work of Zhang and Rundell [36] who present a comparison of global sensitivity approaches to a model of the TCR (T-cell receptor)-activated ERK (extracellular-signal-regulated kinase)/MAPK pathway. Their analysis points to key strengths of a number of available approaches.

### Software tools

A number of software tools carry out sensitivity analysis in the context of MCA, including Jarnac (http://www.sys-bio.org) and Gepasi (http://www.gepasi.org). Dynamic sensitivity analysis can be carried out through BioSens, a toolbox for bio-SPICE (http://www.biospice.org).

## Conclusion

Models developed for systems biology typically depend on many parameters. The construction, analysis, interpretation and application of such models can be greatly facilitated by the use of sensitivity analysis. A strong foundation for the use of local sensitivity analysis has been laid in the work on MCA and the S-system framework of BST. Current developments in the area are novel uses of sensitivity analysis for model calibration and analysis, the expanding use of global sensitivity analysis, and the investigation of models displaying dynamic or stochastic behaviour. These methods will continue to be a key component of our modelling toolbox as the systems biology community takes on more complex systems and more ambitious applications.

## Summary

•

*Sensitivity analysis addresses the manner in which system behaviour depends on parameter values.*•

*Global sensitivity analysis addresses wide regions of parameter space using statistical methods, whereas local sensitivity analysis provides a more complete description of behaviour near a specified operating condition.*•

*Sensitivity analysis is a useful predictive tool in medicine and biotechnology, e.g. in drug development and metabolic engineering.*•

*The construction, analysis and interpretation of models with many (poorly known) parameters, typical of systems biology, is greatly facilitated by sensitivity analysis.*•

*Generalizations of standard sensitivity analysis allow investigation of dynamic or stochastic phenomena.*

- © The Authors Journal compilation © 2008 Biochemical Society