Quantifying robustness of biochemical network models

Abstract

Background

Robustness of mathematical models of biochemical networks is important for validation purposes and can be used as a means of selecting between different competing models. Tools for quantifying parametric robustness are needed.

Results

Two techniques for describing quantitatively the robustness of an oscillatory model were presented and contrasted. Single-parameter bifurcation analysis was used to evaluate the stability robustness of the limit cycle oscillation as well as the frequency and amplitude of oscillations. A tool from control engineering – the structural singular value (SSV) – was used to quantify robust stability of the limit cycle. Using SSV analysis, we find very poor robustness when the model's parameters are allowed to vary.

Conclusion

The results show the usefulness of incorporating SSV analysis to single parameter sensitivity analysis to quantify robustness.

Background

Complex molecular networks mediate intracellular signalling events. These systems must operate reliably under vastly different environmental conditions that can cause large changes in the internal "parameters" of the system. The notion of robustness in biological systems has received considerable interest in the literature recently. By saying that a system is robust we imply that a particular function or characteristic of the system is preserved despite changes in the operating environment of the system. For example, by means of a computer model, Barkai and Leibler demonstrated that the adaptation mechanism found in the chemotactic signalling pathway in Escherichia coli is robust [1]. This was later confirmed experimentally [2]. A model of segment polarity network in Drosophila embryos was also found to be insensitive to variations in kinetic constants that govern its behaviour [3]. A similar approach was later used to show that a core neurogenic network in Drosophila successfully formed three test patterns across a wide range of parameter values [4] leading Meir et al. to propose that the ability to resist parameter fluctuations may be essential for gene network evolutionary flexibility.

Since the signalling pathways are robust, we should expect that mathematical models that attempt to explain these networks also be robust to parameter variations. This has long been appreciated. For example, Savageau, in [5], argues for parameter sensitivities as a means of evaluating the performance of biochemical systems. More recently, Morohashi et al. propose that robustness of a model to parameter variations be used as a criterion for determining plausibility between different models [6].

If we are to use robustness as a means of evaluating the quality of a model, we need objective measures of this robustness. One common technique is through parameter sensitivities. For simple systems, the sensitivity of a model of a network to individual parameters can be evaluated analytically [5, 7]. For more complex networks, it can be determined computationally by repeated simulation varying one parameter while holding all others fixed; [3, 8]. This single parameter sensitivity is also useful for testing robustness of a biochemical network in the laboratory. For example, it is by systematically varying the concentration of the chemotaxis-network proteins in E. coli and determining their effect – or lack thereof – on the precision of adaptation that Alon et al. determined the robustness of this system [2].

Single parameter insensitivity is necessary for a robust network, but may not be sufficient owing to interactions between several parameters. This is particularly true in vivo where many different system parameters will differ from their "nominal" values simultaneously. The tools available for quantifying this multiparametric uncertainty are more limited. Systematic changes of many parameters at a time suffer from an exponential increase in the number of parameters that need to be changed. This "curse of dimensionality" makes varying more than a handful of parameters simultaneously to assess parameter sensitivity impractical. For this reason, sensitivities for several parameters have been traditionally addressed through computer simulations based on Monte Carlo methods [9] – randomly varying all parameter in the model [1, 4]. However, because of their reliance on random methods, Monte Carlo methods cannot guarantee robustness. In this paper we suggest an alternative method, originally developed for use in analysing robust stability in man-made automatic control systems.

The need for robust systems has been one of the primary concerns of control engineering. In fact, one of the earliest motivations for the study of feedback control systems was the need to create robust telephone networks out of the highly variable vacuum tubes of the day. More recently, powerful tools for analysing the robustness of networks have emerged. In this paper we propose that one of these computational tools, known in control theory as the structural singular value (SSV) is of particular interest for biological networks [10]. We do this by contrasting single and multi-parameter sensitivities of a model of an oscillating biochemical network. We describe this model next.

Model of an oscillating biochemical network

In [8], Laub and Loomis propose a model of the molecular network underlying adenosine 3',5'-cyclic monophosphate (cAMP) oscillations observed in fields of chemotactic Dictyostelium discoideum cells. The model, based on the network depicted in Fig. 1, induces the spontaneous oscillations in cAMP observed during the early development of D. discoideum.

Figure 1

Laub and Loomis model. In their model of the aggregation network, pulses of cAMP are produced when adenlylate cyclase (ACA) is activated after the binding of extracellular cAMP to the surface receptor CAR1. When cAMP accumulates internally, it activates the protein kinase PKA. Ligand-bound CAR1 also activates the MAP kinase ERK2. ERK2 is then inactivated by PKA and no longer inhibits the cAMP phosphodiesterase REG A. A protein phosphatase activates REG A such that REG A can hydrolyse internal cAMP. When REG A hydrolyses the internal cAMP, PKA activity is inhibited by its regulatory subunit, and the activities of both ACA and ERK2 go up. Secreted cAMP diffuses between cells before being degraded by the secreted phosphodiesterase PDE.

In this model, changes in the enzymatic activities of these proteins are described by a system of seven non-linear differential equations:

where the state variable x = [x₁,...,x₇] represents the concentrations of the seven proteins: x₁ = [ACA], x₂ = [PKA], x₃ = [ERK2], x₄ = [REG A], x₅ = [Internal cAMP], x₆ = [External cAMP] and x₇ = [CAR1] and the fourteen different k _i represent system parameter values. It was shown numerically in [8] that spontaneous oscillations appear at the nominal parameter values found in Table 1. Note that because there are typographical errors in the original paper, the values being used here for the nominal parameters were obtained directly from the authors of [8].

Table 1 Parameter values For each parameter, k _i denotes the nominal value and the units involved.

This particular model is primarily concerned with describing self-sustaining oscillations in the biochemical system. From Fig. 2, it is clear that at the nominal parameter values of the model, this is achieved. We seek to determine whether the system is robust – that is, if we change these kinetic parameters, will the systems oscillatory behaviour persist? We next present two possible means, based on whether parameters are changed one at a time or in groups.

Figure 2

Oscillations at nominal parameter values. Plot of the concentration of each of the seven variables as a function of time. This figure shows the oscillatory behaviour seen in all the variables.

Methods

Single parameter robustness: Bifurcation analysis

Self-sustained oscillations such as those being modelled here appear as stable limit cycles in trajectory of the underlying dynamical system [11]. The existence and stability of these limit cycles may change under parametric perturbations. Whenever a stable periodic solution loses stability as we vary the underlying parameters of the system and this solution transitions to another qualitative solution – for example, a steady-state equilibrium – we say that the system undergoes a Hopf bifurcation. It is therefore possible to use bifurcation theory as a means of quantifying the robustness of this oscillatory network model [12, 13].

Using the bifurcation analysis package AUTO [14], it is possible to produce one-parameter bifurcation diagrams for each of the model parameters k _i . These diagrams illustrate the steady-state behaviour of the systems as the parameter values are changed. Suppose that Hopf bifurcations occur at k _i and _i so that (stable) limit cycles occur for the range (k _i , _i ). Both the size of this interval as well as the proximity of the nominal parameter value to either boundary are measures of the robustness of the system. To compare the robustness of the system to the different parameters, we can define a series of parametric robustness measures. We define the degree of robustness (DOR) for each parameter k _i as:

It is straightforward to see that this value is always between zero and one. The former indicates extreme parameter sensitivity whereas the latter implies large insensitivity.

Bifurcation diagrams provide an excellent means of determining the robustness of systems to single parameter perturbations. We next describe a method for analysing and quantifying robustness to simultaneous changes in several parameters.

Multiparametric robustness: Structural singular values

As in biological networks, engineering systems must also be robust to variations in the parametric values of its components. Developing tools for the analysis and design of robust automatic control systems has been an area of active research during the last two decades in control theory. One of the most powerful frameworks for measuring robustness known is the structural singular value (SSV) which is due to Doyle and co-workers [15].

We first define and illustrate the use of the SSV to quantify robustness by means of a simple example. Suppose that a system is described by the first-order differential equation

= ax

where the constant parameter a is uncertain, but is assumed to lie in the region a ∈ (a, ).

We would like to know when this system is robustly stable; that is, it is stable for all possible parameters. The differential equation can be rewritten as

= a₀x + b₀v     (1)

where a₀ = ( + a)/2, b₀ = ( - a)/2, v = δx and δ ∈ (-1,1). In Eqn. (1), the term a₀ represents the nominal system description. Clearly, for the system to be stable we need a₀ < 0. The variable δ represents all possible uncertainty in the parameter of the system, whereas b₀ dictates the uncertainty's effect on the nominal model. We would like to maintain system stability no matter what the value of δ happens to be. Note that the system can be redrawn in the form of a feedback loop as in Fig. 3A. From the small gain theorem [16], it is known that if the nominal system is stable (a₀ < 0), the uncertain system will remain so whenever the gain around the loop is less than one. That is, if we denote the system transfer function – the ratio of Fourier transforms of output over input – as

Figure 3

Structured singular value analysis framework. In (A) we show how the uncertain system from Eqn. (1) can be represented in a feedback interconnection involving a nominal system and the uncertainty δ. The signal v = δx provides the feedback. (B) For general systems, uncertainty can also be expressed as a feedback connection of a nominal model and an uncertainty matrix Δ. In this case, the signals u and y provide the interconnections. The transfer function G(ω) = C₀ (iωI - A₀)^-1 B₀ + D₀.

(here i refers to the complex number and ω is the angular frequency) then the system is stable provided that 1 - G(ω)δ ≠ 0, for all frequencies ω. Equivalently, the amount of uncertainty that the system can tolerate is given by

Thus, the function μ = |G(ω)| serves as a (frequency-dependent) measure of the amount of parameter uncertainty that the system can tolerate. In particular, if the size of the uncertainty δ is always less than 1/μ, then the system is robustly stable. In this example, since |δ| < 1, robust stability is guaranteed whenever b₀ < |a₀|.

It is clear that this simple model does not require extensive analytic tools to determine robust stability. Nevertheless, the procedure above can be generalized to systems of the form

where A₀, is a matrix describing the nominal model, B₀, C₀, and D₀ are matrices of appropriate dimensions describing the way that the uncertain parameters affect the nominal model. This uncertainty is modelled by the matrix Δ which is unknown, but is assumed to belong to the set [10]. The signal u = Δy completes the feedback loop as shown in Fig. 3B. For these systems, the appropriate measure of robustness is now given by the structural singular value (SSV)

μ_Δ (G) = (min{||Δ||:Δ ∈ , det(I - G Δ) ≠ 0})^-1     (3)

Here, G(ω) = C₀ (iωI - A₀)^-1 B₀ + D₀ and I is the identity matrix of appropriate dimensions.

Since we are interested in the robustness of the oscillatory property of this system, it is natural to use the SSV to quantify the robust stability of the limit cycle. However, in order to use the SSV tool, the original perturbed system must be transformed into a framework consisting of a nominal linear time invariant (LTI) system interconnected with a perturbation matrix. For the case of an oscillatory non-linear model, this involves several steps, which we outline next.

Determining the limit cycle: Harmonic balance method

The first step is to obtain a mathematical expression for the limit cycle oscillation. The harmonic balance method can be used [17]. The basic idea is to represent the limit cycle by a Fourier series with unknown coefficients (a_n,i, φ_n,i) and period T:

The non-linear differential equation can be used to set up a series of algebraic equations that the coefficients must satisfy. These equations can be solved using numerical packages such as Mathematica or Maple. Depending on the particular form of the limit cycle, a small finite number of coefficients can be used. We can denote this periodic solution as x*(t).

Linearization

The non-linear differential equation must now be linearized about this periodic orbit [17]. Writing the state vector x(t) as

x(t) = x*(t) + x_δ (t)

then the local behaviour of the non-linear system is governed by that of the linearized system:

_δ (t) ≅ J (x*(t))x_≈ (t)

where J is the Jacobian matrix of the vector field f. Note that since the linearization is performed about a periodic orbit, the linear system is periodic.

Restructuring into nominal/uncertainty systems

The Jacobian matrix includes all uncertain parameters. At this point we need to separate the system into a nominal model and a feedback interconnection that involves all parametric uncertainty. We first write each parameter as

where k _i is the nominal value and δ _i is the relative amount of perturbation in the i^th parameter. We now separate the Jacobian matrix as

J (x*(t)) = A₀ (t) + B₀ (t) Δ C₀ (t)     (4)

where A₀(t) is the Jacobian matrix with all parameters at their nominal value, and Δ is a diagonal matrix containing all the uncertainties δ _i . Let y(t) = C₀(t) x_δ(t) and u(t) = Δy(t), the system is now of the form of Eqn. (2) (with x(t) replaced by x_δ(t)).

Discretization

The system can be discretized by sampling the state and output with sampling period h = T/n, where n is a positive integer and assuming that the inputs are piecewise constant; this is also a standard technique in control engineering [18]. In particular, a linear continuous-time system governed by Eqn. (2) gives rise to the discrete-time, linear system

ξ (k + 1) = A _d (k) ξ(k) + B _d (k) v(k)

η (k) = C _d (k) ξ(k)

where A _d (k) = Φ (kh + h,kh), B _d (k) = Φ (kh + h, s) B₀ (s)ds, C _d (k) = C₀ (kh), and Φ (t, τ) is the transition matrix of A₀ (t) [19]. The discretized signals are v(k) = u(kh), η(k) = y(kh), and ξ(k) = x(kh). Periodicity of A _d and B _d is preserved due to the periodicity of the transition matrix. Moreover, it is not difficult to confirm that A _d , B _d and C _d are periodic with period n. The uncertainty matrix after discretization is now Δ _d . The discretization step is illustrated in Fig. 4A.

Figure 4

Discretization using sampling. (A) This diagram illustrates the use of sampling (S) and first-order-hold (H) to discretize a continuous-time system. The sampling circuit's output is equal to the inputs at the sampling times. The first-order-hold circuits generate a piecewise-constant signal equal to the inputs. If we introduce two copies of these circuits into the loop of Fig. 3B, and group the subsystems as shown here, the effect is to generate a discrete-time nominal model G _d as well as a discrete-time uncertainty structure Δ _d . (B) The validity of this approximation will depend on the value of the sampling period, h = T/n, chosen. For the Laub & Loomis model, the error for values of n greater than 8 is negligible. A comparison of the responses of the non-linear (solid blue), linearized continuous-time (dashed green) and discrete-time (dotted red) systems when n = 16 is shown, where we have plotted x₁ as a function of time for all responses. The latter two have been superimposed onto the nominal limit cycle (x*(t)) computed using the harmonic balance method.

Lifting

The final step in preparing the system for SSV analysis is to transform the periodic, linearized system into an equivalent time-invariant one. The technique for this is known as lifting [18]. Rather than giving the general formulae, it is easier to illustrate the general principle with an example.

Suppose that a discrete-time system with state variable ξ, input v, and output η obeys the difference equation

ξ (k + 1) = a(k)ξ(k) + b(k)v(k)

η (k) = c(k)ξ(k)

where the time varying coefficients a(k), b(k) and c(k) are all periodic with period two. Calculating the state variable and output step-by-step leads to:

for any integer p. By defining "lifted" inputs and outputs

and considering the system state only at the even time points (

(p) = ξ(2p)) we arrive at an equivalent time-invariant system.

The lifting technique has been illustrated above for a discrete-time system with period two; however, it can be applied to systems with arbitrary period – though the corresponding formulae are considerably more complicated; see [18].

Computation of SSV

There is considerable literature in control theory on the computation of the SSV; see for example [20–22]. For general classes of uncertainty, computing μ_Δ is known to be NP-hard [21]. Typically, given the feedback loop consisting of G and Δ we compute upper and lower bounds for the SSV [15]. The lower bound is exactly equal to μ_Δ [15]; unfortunately computing this lower bound involves a search over a non-convex set and therefore may converge to local optimums that are not global. In contrast, the upper bound can be rewritten in terms of a convex optimisation problem, so that a global minimum can be obtained. However, this upper bound is, in general not tight. A software package is commercially available that can compute μ^upperand uses a power algorithm to attempt to compute μ^lower[22].

Results

Single parameter robustness

The robustness of Laub and Loomis's oscillatory model was first analysed by means of single-parameter bifurcation diagrams. Four typical diagrams are shown in Fig. 5. The activity of internal cAMP (x₅) is plotted as a function of the variation of each parameter. We use internal cAMP in the diagram as it is the element that is usually observed experimentally [23]. In each of the diagrams, there are three types of solutions: stable steady state, unstable steady state and limit cycle oscillations.

Figure 5

Single parameter bifurcation diagrams. In these diagrams, obtained using AUTO [14], we plot the steady state activity of internal cAMP (x₅) as a function of four individual parameters (k₁, k₂, k₈ and k₁₀). In the diagrams, a stable steady state is represented by solid line; an unstable steady state is represented by dashed line; and a stable limit cycle is represented by a pair of solid circles with the upper indicating the maximum value of amplitude and the lower indicating the minimum value of amplitude. The transition from stable steady state to stable limit cycle or vice versa is called a Hopf bifurcation, which is indicated by a solid square in the plot. This type of bifurcation is caused by the appearance of a pair of pure imaginary eigenvalues of the Jacobian matrix of non-linear system. A star on the axis indicates the nominal value of each parameter.

These diagrams illustrate that Hopf bifurcations occur for each parameter; that is, the oscillatory behaviour exists only in a limited range of parameters around the nominal value. For each of these parameters, the respective intervals and values for degree-of-robustness are found in Table 1.

Structural singular value

From the numerical simulation (Fig. 2) of the non-linear model, we observed that the oscillatory curves did not deviate greatly from a simple harmonic oscillator plus a constant offset. Thus, to obtain an analytic expression for the periodic orbits we assume that the state variables are expanded into Fourier series containing only the fundamental and constant terms:

for each of the seven states. Since it is the relative phase shift between each state variable that is relevant, we assume that θ₁ = 0. The substitution of the Fourier series into the original equations leads to a series of real algebraic equations for the coefficients (not shown) whose solution was obtained using Mathematica. This leads us to obtain the corresponding periodic solutions where the values of A_0,i, A_1,iand θ _i are found in Table 2. The period T is approximately 7.31 minutes. This analytic solution matches well with the numerical simulation except for an arbitrary phase shift, which does not affect the shape and location of the limit cycle in the phase space and can thereby be ignored (not shown).

Table 2 Fourier coefficients of the nominal periodic solution Values obtained for our Fourier series opproximation of the oscillatory behaviour.

Following our prescribed methods, we next linearized the system. The Jacobian matrix is obtained and was decomposed as in Eqn. (4) to obtain:

The matrix B₀ (t) = {B_i,j} where

Similarly, the matrix C₀ (t) = {C_i,j} where all coefficients are zero except for the following:

Finally, D₀ = 0 and the perturbation structure is given by

Δ = diag {δ₁,δ₂,δ₂,δ₃,δ₄,δ₅,δ₆,δ₆,δ₈,δ₈,δ₉,δ₁₀,δ₁₀,δ₁₁,δ₁₂,δ₁₃,δ₁₄}

Note that since the nominal trajectory is periodic, the matrix functions in the nominal description are also periodic. Note also that in the uncertainty matrix, Δ, some uncertainties are repeated (δ₂,δ₆,δ₈ and δ₁₀) while δ₇ is missing.

The system was then discretized and lifted following the procedure outlined above. A comparison of the system response for each of the approximations at the nominal parameter values is given in Fig. 4B.

For the sampling time we tried various values of n but found negligible differences for values above eight. Finally, we computed the bounds on the SSV. Once again, we found that the values of these two bounds were not affected much by the sampling frequency provided that n is greater than eight. The upper bound was successfully computed using [22]. The maximum over all frequencies is approximately 12.06. However, the high dimension of system causes convergence problem during the computation of μ^lowerusing this package. To obtain an acceptable lower bound, we calculate the spectral radius at each frequency. This gives us a lower bound for μ^lower. The plot of the bounds for the SSV when n = 16 is shown in Fig. 6. We can use μ^lowerto obtain a conservative region for robust stability. The highest value over all frequencies for μ^loweris approximately 2.636.

Figure 6

Structured singular value bounds. Shown are the lower (dashed) and upper (solid) bounds for the structural singular value as a function of the frequency. Since we are interested in the largest value of μ over all frequencies, we find that the maximum value of μ is between 2.636 and 12.06.

Discussion

Recent years has seen an appreciation that key cellular properties are robust to variations in individual parameter values. Based on the topology of many of these networks, this should not be surprising. Feedback – both negative and positive – control systems are ubiquitous in most biological networks [24] and one of the reasons for using feedback is that it reduces sensitivity of a system's behaviour to its parameter values.

In modelling biological networks, it is important that this robustness also be in evidence. The particular behaviour being characterized by the model should not rely on precise values of the model's parameters – for example, reaction rate constants or protein concentrations. In particular, a precise measurement of these constants is difficult whereas protein concentrations will vary from one cell to another or throughout the lifetime of any individual. Deviations from the nominal model parameter values should not result in a loss of the network's performance; thus, parameter sensitivity can be used to validate mathematical models of biochemical system. That is, the more insensitive the system response is to the accuracy of the parameter, the more faith we should have in the model [25].

In looking at certain classes of behaviour, where qualitative changes in the stability of the system are possible, bifurcation diagrams provide an elegant means of evaluating robustness. For example, in evaluating the robustness of the model of Laub and Loomis, of primary importance is determining whether the oscillatory behaviour will persist if the parameter values are altered. This qualitative difference in performance – from limit cycle oscillations to constant steady states – can be quantified and compared across parameters or from one model to another. Once the robustness of the oscillatory behaviour is established, further investigations of the robustness of some of the oscillatory features, for example frequency and amplitude can further be evaluated.

From the bifurcation diagrams obtained for each of the fourteen parameters, we know that oscillations exist only in a limited range around the nominal value. We find the system to be quite sensitive to variations in k₂, k₄, k₁₀ and k₁₄ and mostly insensitive to the others. Single-parameter bifurcation analysis also shows that the amplitude of the oscillation is greatly affected by the variation of 9 parameters (k₁, k₂, k₄, k₆, k₇, k₁₀, k₁₁, k₁₂, and k₁₄).

Based on the SSV stability to interpret multiple parameter sensitivity, we can conclude that robust stability of the periodic orbits will be maintained, provided that

Since the uncertainty matrix consists only of diagonal entries, this bound applies to each of the individual parameters. Thus, we can guarantee that the system will be robustly stable provided that no single parameter differs more than 8.3% from its nominal value.

In our analysis we found a large gap between μ^upperand our lower bound for μ. As we later show, for this system the upper bound is fairly tight, as we are able to obtain a destabilizing perturbation of size 9%. For general biological models, a robustness measure based on the upper bound μ^uppermay also be more appropriate. Robustness bounds for systems in which arbitrarily slowly-time-varying parameter values are allowed are known [26]. For these systems it has been shown that the bounds converge as the time-variations approach zero to the upper bound μ^upper[26]. Since many of the parameters in models of biochemical networks represent features that will vary over time, such as enzyme concentrations, this number may therefore be more indicative of the model's true robustness.

The ability to consider the effect of time-variations on the robustness of the system is one great advantage of the SSV over other methodologies. One drawback of the SSV approach compared to the bifurcation theory is that it does not provide the precise combination of parameters that destabilizes the system – only its size. Also, since the upper bound is only sufficient to guarantee robustness, this number may, in general, give an overly conservative notion of robustness.

It must be emphasized that the SSV approach denoted here is based on the linearized model of the system. For some classes of systems this linearization may not be possible – in this case, the linear SSV approach documented here will not be applicable. However, for most models used to describe biochemical reactions, this should not be a problem.

Because we are concentrating exclusively on the local stability of the linearized model, important parameters of the oscillatory behaviour such as robustness of the frequency and amplitudes of oscillation are not evaluated. Also, the effect of parameter variations on the equilibrium orbit are omitted. In particular, varying the kinetic parameters will change the behaviour of the system in two different ways: the equilibrium periodic orbit will change and the stability of deviations about this orbit will also change. The SSV allows one to quantify the robustness of the second of these two effects. It does not say anything directly regarding the effect of parameter variations on the equilibrium periodic orbit. One way of bounding the effect of these parameter changes is to write the original differential equation as

(t) = f(x,k)

where k = k₀ + δ is the set of kinetic parameters with nominal values k₀. If the nominal periodic orbit (when δ = 0) is given by

x*(t) = x(t) - x_δ (t)

then, linearizing about this orbit yields

_δ(t) = A(t)x_δ (t) + v

where

δ is a constant vector that includes the effect of this parametric uncertainty. Thus, the system can be considered as being perturbed by a constant input signal v. Provided that the homogeneous system is exponentially stable (and this is guaranteed by the existence of a stable limit cycle) and that v is not "too large", the perturbed system's state will remain in a neighbourhood of the origin if the f(x,p) in the original equation is reasonably well behaved in k. Detailed bounds and conditions on f are given in Theorem 5.1 of [17], though it should be emphasized that these bounds tend to be overly conservative in practice.

To illustrate the local nature of the SSV analysis for this system, we perturbed the system parameters by varying amounts. The particular parameters were either increased or decreased so as to bring them closer to the Hopf bifurcation. For example, the nominal value of k₁ is closer to k₁ than to ₁ so that we reduced k₁ whereas k₄ is closer to ₄ than to k₄ so we increased k₄. In Fig. 7 we show the response of these systems to changes of 7% and 9% both for the linearized system – where the linearized response has been superimposed on the nominal limit cycle (Fig. 7A) and the original non-linear system (Fig. 7B). For the smaller value, the linearized response is stable and we see that, after a transient, the response settles to the nominal limit cycle. We also see this same behaviour in the response of the non-linear system with this level of parameter perturbation. For a 9% change in the parameters, however, the linearized system is unstable. We see this as a deviation from the nominal limit cycle. In the non-linear system's response, this translates into an end to the stable limit cycle. The response does not "blow up" but instead settles into a fixed point.

Figure 7

Comparison of linearized and non-linear responses for varying parameters. We show the response of the system to two levels of perturbations of size 7% and 9%. (A) The linearized perturbed system superimposed on the approximate nominal limit cycle obtained using the harmonic balance method. (B) The response of the non-linear system. For the perturbation of size 7% the linearized system is stable and so we see the responses settle into a stable limit cycle. For 9%, the linearized system is no longer stable, so that the response of the linear system increases unboundedly, as the trajectory moves away from the nominal equilibrium. In the non-linear response we see this manifested as a transition to a stable steady state.

This example illustrates how a robustness analysis of the linearized system can be used to deduce the robustness of the original non-linear system, as it shows that when the linearized system is unstable, the desired behaviour of the non-linear system will no longer be present. This example also points out the fact that the upper bound μ^upper≅ 8.3% is not overly conservative for this system as we were able to produce a destabilizing perturbation of size 9%.

Finally, we note that multiparametric robustness analysis considered here is based on local properties of the dynamical system, since we are evaluating the robustness of the linearized model. Extensions to the non-linear model are the subject of active investigation [27].

However, it is by combining the robustness analysis of both single and multiple parameters, we can obtain a more thorough understanding of the region of stability of the periodic solution in the high dimensional parameter space and use this to improve upon the model. In this particular example, we find that the system's robustness is governed by several "robustness limiting" parameters, k₂, k₄, k₁₀ and k₁₄.

Conclusions

Determining the robustness of mathematical models of biological systems is important for several reasons. First, there is growing evidence that many aspects of the networks being modelled have evolved in such a way so that they are robust as this allows them to tolerate natural variations in the environment. Thus, faithful models should replicate this robustness. Second, robustness of the models provide a means of validating model quality since the performance of the models should not rely on precisely tuned parameter values that are impossible – or at best – difficult to measure exactly.

In this paper, we illustrated the use of two tools developed in dynamical systems theory and control engineering to assess robustness quantitatively. For an example, we considered an oscillatory molecular network model due to Laub and Loomis that aims at describing oscillatory behaviour in cAMP signalling observed in the social amoeba D. discoideum. This behaviour appears as a stable limit cycle of the equations describing the model. We have evaluated the degree to which this limit cycle is robust to variations in all the system parameters.

The robustness of the oscillatory behaviour to single parameter variations was quantified using bifurcation analysis. Using the bifurcation analysis software tool AUTO we determined that single parameter changes as small as 20% from the nominal value can cause the limit cycle to disappear and a stable equilibrium to appear. In addition to the stability robustness, AUTO is also able to evaluate the sensitivity of the amplitude of the oscillation to these parameter changes.

To investigate the robustness of the model to simultaneous changes in parameter values, the structured singular value (SSV) analysis tool was used. Once the system was in the correct framework for SSV analysis, we were able to determine that the system can only tolerate very small changes in the parameter values – in the order of 8% – if we allow these parameters to vary with time arbitrarily slowly.

Finally, it is important to note that to understand completely the robustness properties of a model, it is appropriate to combine single and multiple parameter sensitivity analyses.