- Email: [email protected]
Sensitivity analysis of oscillatory (bio)chemical systems
Computers and Chemical Engineering 29 (2005) 663–673
Sensitivity analysis of oscillatory (bio)chemical systems Daniel E. Zaka , J¨org Stellingb , Francis J. Doyle IIIc,∗ b
a Department of Chemical Engineering, University of Delaware, Newark, DE 19716, USA Max Planck Institute for Dynamics of Complex Technical Systems, D-39106 Magdeburg, Germany c Department of Chemical Engineering, University of California, Santa Barbara, CA 93106, USA
Received 22 September 2003; received in revised form 8 January 2004 Available online 7 October 2004
Abstract From the cell cycle to circadian rhythms, oscillatory processes are fundamental to biology. Emerging from nonlinear dynamical interactions, oscillatory mechanisms are best understood through mathematical modeling. Ordinary differential equations (ODEs) are one framework in which the complex interactions giving rise to biological oscillations may be modeled. Key to ODE models are the model parameters that determine whether or not oscillations will occur, and the period and amplitude of the oscillations when they do. Sensitivity analysis is a means to acquire insight about the importance of the model parameters. Sensitivity analysis of oscillatory systems provides unique challenges and must be addressed carefully. In the present study, we describe a method for determining the sensitivity of the period to the model parameters that is straightforward to implement and interpret. We apply this method to a model for circadian rhythms, and obtain results suggesting a link between network structure and parameter sensitivity. © 2004 Elsevier Ltd. All rights reserved. Keywords: Circadian rhythm; Mathematical modeling; Sensitivity analysis
1. Introduction Rhythmic processes are fundamental to biology, as is revealed by literature searches on two key oscillatory biological processes, ‘circadian rhythm’ and ‘cell cycle’, that yielded 43,057 and 61,556 articles, respectively (http://www.ncbi. nlm.nih.gov: December 22, 2003). Circadian rhythms are variations in activity (physiological and cellular) with nearly twenty-four hour periods that, for the fly and the mouse, arise from cellular genetic networks containing delayed feedback mechanisms (Hastings, 2000). The cell cycle is the process by which a single cell divides to become two, and is comprised of a complex dynamic nonlinear network of protein interactions (Tyson, Csikasz-Nagy, & Novak, 2002). The oscillations in these and other cellular processes arise from biochemical reactions, and thus general principles that have been developed for chemical oscillators apply to biological oscillators. One ∗
Corresponding author. E-mail address: [email protected] (F.J. Doyle).
0098-1354/$ – see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2004.08.021
principle is that computational approaches are necessary to fully understand and explain their behavior (Rabitz & Edelson, 1985; Goldbeter, 1996; Goldbeter, 2002). The approaches employed for the mathematical analysis of biological oscillators have paralleled those used for chemical oscillators, where both ordinary differential equation (ODE) formulations, derived from mass action kinetics, and discrete stochastic formulations (Gillespie, 1976) have been employed, with the ODE approaches being predominant (Rabitz & Edelson, 1985). Interest in stochastic models of biological oscillators (Barkai & Leibler, 2000; Gonze, Halloy, & Goldbeter, 2002b; Vilar, Kueh, Barkai, & Leibler, 2002; Zak, Doyle, Vlachos, & Schwaber, 2001), and of biological systems in general, is growing, however, and has been strengthened by the publication of several elegant experimental studies that demonstrated the stochastic nature of some biochemical reactions at the single-cell level (Levsky, Shenoy, Pezo, & Singer, 2002; Elowitz, Levine, Siggia, & Swain, 2002), reviewed in (Rao, Wolf, & Arkin, 2002). In spite of this, new models of biological oscillators in the ODE framework continue to be developed (Leloup & Goldbeter,
664
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
Nomenclature A B f p p0 S(t) Sτ S sτ S c (t) sij sτj sτsj v vp x
ns × ns state Jacobian matrix ns × np parameter Jacobian matrix vector of state derivatives of size ns vector of parameters of size np nominal parameter values ns × np sensitivity matrix period sensitivity vector scaled period sensitivity vector ns × np cleaned out sensitivity matrix sensitivity of the ith state to the jth parameter sensitivity of the period to the jth parameter scaled sensitivity of the period to the jth parameter input direction from SVD perturbation direction vector state vector of size ns
Greek letters ps scaled perturbation strength τ s scaled period deviation (oscillator precision) ϕ perturbation scaling factor σ singular value τ period of the oscillator
mechanistic insight and to design informative experiments. Sensitivity analysis is one technique to investigate the importance of parameters. Sensitivity analysis has been applied in a few cases to the analysis of biological systems, with the objectives of mechanism discrimination on the basis of model sensitivity (Morohashi et al., 2002; Savageau, 1971), experimental design and parameter estimation (Leif & Jorgensen, 2001; Schlosser, 1994), and the relationship between sensitivity and identifiability (Stelling & Gilles, 2001). In the present work, we consider the parametric sensitivity analysis of ODE models of biological oscillators, where the parameters determine both the existence of and characteristics of the oscillations. Since the period of biological oscillators is a key aspect of their physiological significance (for example, the time-keeping nature of circadian rhythms), we specifically consider the parametric sensitivity of the period. Direct application of sensitivity analysis to oscillating systems gives rise to secular terms (Larter, 1983), and thus must be done with care. In the present work, we first discuss the standard methods for the sensitivity analysis of oscillatory systems, and then present a novel method based on singular value decomposition (SVD) that is easier to interpret and implement than the common methods. We conclude with a case study of an ODE model for circadian rhythms. 2. Sensitivity analysis of oscillatory systems
2003; Forger & Peskin, 2003), and stochastic and deterministic simulations are often used as complementary, rather than competing, methods by individual groups (Gonze, Halloy, & Goldbeter, 2002a; Gonze, Halloy, Leloup, & Goldbeter, 2003; Leloup & Goldbeter, 1999; Novak, Pataki, Ciliberto, & Tyson, 2001; Sveiczer, Tyson, & Novak, 2001; Vilar et al., 2002). Generally, the ODE framework is used because of its ease of simulation and analysis with available software tools (Maly & Petzold, 1996; Shampine & Reichelt, 1997), while the stochastic framework is used to explore fluctuations. Interestingly, there have been a number of studies of biological oscillators where predictions made using the ODE framework have held up in a stochastic framework, with the effect of stochastic noise largely being the widening of the limit cycle (Gonze, Halloy, & Goldbeter, 2002a; Gonze et al., 2003). Thus, while experimental results have demonstrated that the stochastic nature of some biological reactions cannot be neglected, ODEs, and tools to analyze them, are likely to continue to be of importance to the computational biology community. Given that oscillating chemical systems typically consist of larger numbers of molecules and potentially faster reactions, we also expect ODEs to continue to play an important role in the study of chemical oscillators. ODE models of biological processes often involve many parameters, and the importance of these parameters in determining system behavior must be assessed in order to gain
In the present section, we first describe basic principles of sensitivity analysis that are applicable to any ODE system. This is followed by a presentation of aspects of sensitivity analysis that are specific to oscillatory systems. 2.1. General principles of sensitivity analysis Consider the ODE system: x˙ = f (x, p)
(1)
where x(t) is a vector of ns states, p is a vector of np model parameters, and f is a column vector of the state time derivatives. Assuming that a solution to Eq. (1) exists, the sensitivity matrix of the system, S(t), that describes how variations in the parameters near the local point in parameter space, p0 , influence the state trajectories, may be defined: ∂x (2) S(t) ≡ ∂p x=x(t,p0 ),p=p0 where S(t) is composed of individual sensitivities of each state to each parameter (sij ). The simplest way to calculate S(t) is by finite differences. For a single parameter pj , x(t, pj + pj ) − x(t, pj ) ∂x ≈ (3) ∂pj x=x(t,p ) pj 0
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
While Eq. (3) may be used to estimate the sensitivities, it is computationally tedious and frequently inaccurate, given that it may lead to numerical instabilities (Rabitz & Edelson, 1985). An alternative approach is to differentiate Eq. (1) with respect to the parameter vector p, giving (Khalil, 1992): S˙ = A(t, p0 )S + B(t, p0 )
2.2. Oscillatory systems The present work is concerned with systems that are periodic in time, for which: x(t + τ) = x(t)
(6)
where τ is the period of the oscillation. The objective is to determine how τ depends on the model parameters, p, specifically to determine the np × 1 vector of period sensitivities, Sτ : T ∂τ ∂τ ∂τ Sτ ≡ (7) , ,···, ∂p1 ∂p2 ∂pnp S τ is composed of individual parameter sensitivities, sτj , and is a constant vector. Given that x(t) is periodic, it is possible to express each of the states x1 , x2 , . . . xns as Fourier series (Larter, 1983):
n=0
ani cos
2nπt 2nπt + bni sin τ τ
(8)
where ani and bni are Fourier coefficients and are functions of parameters. The time derivative x˙ i can also be calculated: ∞ 2π 2nπt 2nπt −nani sin + nbni cos τ τ τ n=0
n=0
∂ani 2nπt ∂bi 2nπt cos + n sin ∂pj τ ∂pj τ
(10)
Using Eqs. 8 and 9, Eq. (10) can be rearranged to: (5)
The sensitivity matrix S(t) as defined by Eq. (4) is a linear time-varying (LTV) system. By symbolically calculating the state and parameter Jacobians (A and B, respectively) it is possible to integrate x(t) and S(t) simultaneously, given the nominal parameter values, p0 , and the initial conditions for the states, x0 . Unless the initial conditions of the states depend on the parameters, the initial conditions for the sensitivities are zero. This method is known as the direct method (Rabitz & Edelson, 1985), and involves the integration of a system of dimension ns × (np + 1). When np > ns , Eq. (5) may be more efficiently integrated through the use of a Green’s function technique (Rabitz & Edelson, 1985).
∞
∞ n=0
∂f ∂x x=x(t,p0 ),p=p0 ∂f B(t, p0 ) = ∂p x=x(t,p0 ),p=p0
x˙ i (t) =
∞ 2πt ∂τ 2nπt 2nπt i i sij = 2 nan sin − nbn cos τ τ τ ∂pj
(4)
A(t, p0 ) =
xi (t) =
Finally, the sensitivity sij may be calculated, keeping in mind that τ depends on p:
+
where
665
(9)
t sij = − sτj fi + τ
∂xi ∂pj
(11) τ
or in vector notation: t S = − f S Tτ + S c τ
(12)
where S c (t) is the ns × np cleaned out sensitivity matrix (Tomovi´c & Vukobratovi´c, 1972), composed of (∂xi /∂pj )τ , and is periodic in time, as can be seen in Eq. (10). S c (t) contains information about how variations in the parameters affect the shape of the trajectory (Larter, 1983). The above Fourier series-based derivation was given by Larter (1983), but Eq. (11) was originally derived by Tomovi´c and Vukobratovi´c (1972). Eq. (12) clearly demonstrates that, because of the secular term, the raw sensitivity coefficients sij are unbounded in time for parameters with sτj = 0. A similar observation has been made in the metabolic control analysis of autonomous oscillators (Kholodenko, Demin, & Westerhoff, 1997). The secular term arises because the computation of the sensitivities involves a nonuniformly valid expansion of a periodic system (Larter, 1983). Through the Linstedt–Poincar´e method, Larter (1983) has shown that the periodic cleaned out sensitivity matrix, S c (t) is actually the first correction term in a uniformly valid expansion of the system around the nominal parameter set (Larter, 1983). To determine either how variation in the parameters influences the oscillator period (S τ ) or the shape of its trajectory (S c (t)), the raw sensitivity matrix (S) must be broken into its component terms. This is generally accomplished by calculating S τ by another means, and then applying a correction to S(t) (Larter, 1983). Two methods for computing S τ have been described in the literature, and these will be considered in the following section. 2.3. Current methods for determining period sensitivities Edelson and Thomas (1981) have described one approach for calculating S τ . Differentiating Eq. (6) for a particular state, xi , with respect to a particular parameter,
666
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
pj , gives: dxi (t, p) dxi (t + τ, p) = dpj dpj ∂xi ∂t ∂xi ∂xi ∂(t + τ) ∂xi + = + ∂t ∂pj ∂pj t ∂(t + τ) ∂pj ∂pj t+τ ∂xi ∂xi ∂τ ∂xi = + ∂pj t ∂t ∂pj t ∂pj t+τ
Eq. (16) may be simplified by combining with Eq. (11). Firstly: sτj ∂xi (17) txi − xi dt + sij dt = − τ0 ∂pj τ (13)
Rearranging gives:
=
Sij (t) − Sij (t + τ) fi (t)
(14)
In practice, Eq. (14) additionally will contain a transient term that decays over time to yield the correct sτj , even if the initial conditions lie on the limit cycle (Rabitz & Edelson, 1985). This transient results from the fact that these initial conditions are actually dependent on the model parameters (Rabitz & Edelson, 1985). It must be pointed out that computation of sτj by Eq. (14) involves division by the time derivative of xi , fi , that must be zero at some points along the trajectory because the system is oscillatory. For this reason, period sensitivities computed using this method are not defined at all points along the cycle. As described by Edelson and Thomas (1981) themselves, care also must be taken in choosing which state is used to compute the period sensitivities, to avoid both division by small time derivatives and inaccuracies in the difference in the numerator. For the Oregonator system they considered, two out of the three states of the dynamic system proved to be undesirable for computing period sensitivities (Edelson & Thomas, 1981). An alternative integral approach to computing period sensitivities that removes the need to divide by time derivatives (but not does not resolve the issue of initial transients) has been given by Larter (1983): t1 +τ sτj = −
t1
sij dt −
t2 +τ t2
sij dt
(15)
xi (t1 ) − xi (t2 )
Eq. (15) will not exhibit instabilities as long as t1 and t2 are chosen so that xi (t1 ) − xi (t2 ) is not small. This improvement in numerical performance is traded off against an increased demand on the accuracy of the period (τ) used in the integrations, however. This aspect of Eq. (15) arises through the numerator, and is explained in the following derivation. Suppose that the period, τ, used for the integrations in Eq. (15) is actually the true period of the oscillator, τ0 , plus some small error . This gives:
tk +τ
tk
sij dt =
tk +τ0
tk
sij dt +
tk tk +τ0 tk
(∂xi /∂pj )t − (∂xi /∂pj )t+τ = (∂xi /∂t)t
sτj
Since xi and (∂xi /∂pj )τ are both periodic, the following are true for any time tk : tk +τ0 xi dt = Cx
tk +τ0 + tk +τ0
sij dt
(16)
∂xi ∂pj
(18) τ
dt = Cτ
where both Cx and Cτ are constants. The following are reasonable approximations when is small: tk +τ0 + xi dt ≈ Cx tk tk +τ0 + ∂x i
tk
∂pj
τ
dt ≈ Cτ
(19)
xi (tk ≤ t ≤ tk + ) ≈ fi (tk )t + b(tk ) where fi (tk ) and b(tk ) are from the local linearization of xi in time and are periodic. Combining Eq. (16) with Eqs. 17–19, and neglecting all terms of order and higher, gives: tk +τ sτj sτj sij dt ≈ −sτj xi (tk )+ Cx + Cτ − fi (tk )tk (20) τ τ0 0 tk Substituting Eq. (20) into Eq. (15) gives: sτcalc = sτj j
(xi (t1 ) − xi (t2 )) + (/τ0 )(fi (t1 )t1 − fi (t2 )t2 ) (21) xi (t1 ) − xi (t2 )
In general, fi (t1 ) and fi (t2 ) will differ, and thus the error will grow unbounded with increasing t1 and t2 . The rate at which the error grows can be slowed by choosing t1 and t2 so that fi (t1 ) ≈ fi (t2 ). To summarize, the current methods for computing S τ involve calculating sτj for each parameter individually, and require careful selection of the states and times at which the calculations are performed. For the method of Edelson & Thomas (1981), t must be chosen so that fi is not close to zero. For the method of Larter (1983), the integrated sensitivities for at least one state xi (t) must be computed, and t1 and t2 must be chosen to maximize the difference between xi (t1 ) and xi (t2 ), while minimizing the difference between fi (t1 ) and fi (t2 ). A method that can determine all of the period sensitivities at once and without concerns of numerical stability is desirable and is the focus of the following section. 3. Singular value decomposition approach to period sensitivity The novel method for determining period sensitivities that is proposed in the present work is based on the singular value decomposition (SVD) of the raw sensitivity matrix S(t).
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
667
SVD has seen widespread application in multivariable feedback control, to identify strong and weak input directions, as well as to quantify interactions between inputs and outputs (Skogestad & Postlethwaite, 2000). SVD has also been applied to linear and LTV discrete time systems to gain a better understanding of non-minimum phase systems in the time domain (Wang & Lee, 1996) and the reachable subspace of system outputs (Wang & Zachery, 1996). The following derivation is based on the observation that, at large times, the first term on the right-hand side of Eq. (12) will dominate the periodic term S c , giving:
Postlethwaite, 2000), the single non-zero eigenvalue of S τ S Tτ is simply S Tτ S τ . This gives:
t S = − f S Tτ τ
v = αS τ
(22)
The value of t that constitutes sufficiently large will depend on the characteristics of the system being analyzed, specifically the magnitudes of τ, f , and S τ . To be valid, Eq. (22) also requires f = 0, and S τ = 0. Given that x(t) is periodic, the first condition obviously cannot be guaranteed for all points along the trajectory for all components of f . However, the periodicity of the system also guarantees that all components of f will never be zero simultaneously, because that would imply steady-state behavior. Except for some pathological systems (Larter, Rabitz, & Kramer, 1984), τ will be sensitive to at least one parameter, and it is a valid assumption that S τ = 0. The rare cases where it is not will be obvious from inspection of S(t), which will be purely periodic. The SVD of S(t) gives: S = UΣV T
(23)
where Σ is an ns × np diagonal matrix of non-negative singular values, σk , that are the square roots of the eigenvalues of S T S or SS T . V is an np × np matrix whose columns are unit eigenvectors of S T S. U is an ns × ns matrix whose columns are unit eigenvectors of SS T . The vectors which comprise U and V are termed the input and output vectors of S given that: Svk = σk uk
(24)
Thus, given an input along direction vk , the output will be in the direction uk with gain σk (Skogestad & Postlethwaite, 2000). Given that SVD of S(t) requires operations on the matrix S T S, it is convenient to simplify it using Eq. (22): 2 t2 T T 2t (25) S f f S = φ S τ S Tτ τ τ τ2 τ2
where φ(t) = f T f = 0 (since φ ≡ 0 would imply steadystate). It follows that:
STS =
eig(S T S) = φ2
t2 eig(S τ S Tτ ) τ2
(26)
where the expression eig indicates eigenvalues. Since, for any matrix X, the non-zero eigenvalues of XT X are equivalent to the eigenvalues of XXT (Skogestad &
t2 T (27) S Sτ τ2 τ To find the eigenvector v that corresponds to λ, and therefore the most sensitive input direction, the following condition must be satisfied:
eig(S T S) = λ = φ2
S T Sv = λv
(28)
The solution is: (29)
where α is an arbitrary scalar. The solution is readily verified: S T SαS τ = αφ2
t2 S τ S Tτ S τ = αλS τ = λv τ2
(30)
Thus, singular value decomposition of the raw sensitivity matrix at large t will yield a scalar multiple of the vector of period sensitivities, S τ , as the input vector corresponding to the only nonzero singular value. The determination of the magnitudes of the period sensitivities from the SVD now follows. Since eigenvectors are usually normalized to be unit vectors, the scalar multiple α will be:
−1 T α=± Sτ Sτ (31) Thus the input vector corresponding to the largest singular value will be: v = ±
Sτ
(32)
S Tτ S τ
where the relative, but not absolute, signs of the period sensitivities are known. The magnitudes of the period sensitivities can be found using the singular value. Given that:
√ t σ= λ=φ S Tτ S τ (33) τ it follows that: στ Sτ = ± v φt
(34)
In conclusion, the magnitudes and relative signs of the period sensitivities of an oscillatory system may be determined using the SVD of the raw sensitivity matrix at large t, when f = 0, and when τ is sensitive to at least one parameter. This result has an interesting interpretation in terms of the classical application of SVD in that the vector of period sensitivities corresponds to the most sensitive input direction of the raw sensitivity matrix. What constitutes sufficiently large for t will depend on the specific system analyzed, and thus should be determined by considering the convergence of the calculated values of S τ . Also, when Eq. (22) is strictly true, there will only be one non-zero singular value. Eq. (22) is an approximation to a system truly described by Eq. (12), however, and
668
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
Fig. 1. Raw sensitivities of the Goldbeter (1995) model, parameter set A, displaying unbounded growth linear in time as predicted by Eq. (12).
thus there will be several non-zero singular values. At large t, the first singular value, σ1 will dominate the others by several orders of magnitude, and under these conditions σ1 and v1 may be substituted for σ and v in the above expressions. The absolute sign of S τ also must be calibrated, which may be accomplished by a small perturbation to one of the parameters and determining if it increases or decreases the period. Overall, the present method has the advantages over the previous methods in that less care and human intervention is required in selecting states and times for the computations, and that the period sensitivity estimates are guaranteed to converge to their true values. These advantages come at the cost of increased computational requirements, however. 4. Case study: period sensitivity of Goldbeter (1995) circadian rhythm model The present section is concerned with the computation of period sensitivities from the raw sensitivities for the Goldbeter (1995) model for intracellular circadian rhythms.
As mentioned above, circadian rhythms are periodic daily variations of activity that arise from intracellular feedback mechanisms (Hastings, 2000). The Goldbeter (1995) model is a 5-state system (ns = 5) with 18 parameters (np = 18). Two sets of parameter values have been given, with set A being from the original paper (Goldbeter, 1995), yielding a period of 23.7 h. Parameter set B is from a recent publication (Gonze et al., 2002b) and yields a period of 22.2 h. The model equations and parameters are provided in the appendix. Oscillations arise in the system from the general principle of delayed negative feedback. The raw sensitivities for the system were calculated using Eq. (4) by simultaneous integration of the original model and the sensitivity equations (the direct method) in MATLAB (The Mathworks, Inc.) using the function ode15s (Shampine & Reichelt, 1997). The sensitivities were unbounded in time, as predicted by Eq. (12), shown for parameter set A in Fig. 1. The method of Larter (1983) (Eq. (15)) was first used to compute S τ for the Goldbeter (1995) model. In the present calculations, t1 was varied continuously and t2 was chosen to maximize the difference |xi (t1 ) − xi (t2 )| within one τ of time, (t2 < t1 + τ), ensuring that the denominator in Eq. (15) did not approach zero. Although Eq. (15) predicts that the period sensitivities should be constant over time (after the decay of initial transients), they were not, as shown in Fig. 2 for parameter set A. For the parameter vs , for example (Fig. 2a), the period sensitivity showed damped oscillations, eventually settling down to a constant value. This parameter was the exception, as for the remainder of the parameters for parameter set A, and for all of the parameters for parameter set B, the computed period sensitivities grew unbounded in time. This is shown for parameter k2 in Fig. 2b. This unbounded growth was due to the amplification of small errors in the value of τ that was used in the integrations, as discussed above. To obtain actual values of S τ , median values were computed from the time courses. Given values for S τ and Eq. (12), it was possible to calculate the cleaned out sensitivities S c (t). A representative plot
Fig. 2. Period sensitivities to (a) vs (∂τ/∂vs ) and (b) k2 (∂τ/∂k2 ) calculated using the method of Larter (1983).
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
669
Fig. 3. Raw and cleaned out sensitivities (∂MP /∂ks ) for 5-state Goldbeter model, parameter set A. Note that the cleaned out sensitivity is periodic in time as expected.
of a periodic cleaned out sensitivity with its corresponding unbounded raw sensitivity is shown in Fig. 3. The SVD approach also was used to compute the S τ for the Goldbeter model. As predicted by the derivations, at large t, the first singular value, σ1 , dominated the others by several orders of magnitude, demonstrating that the assumption made in deriving Eq. (22) was valid at large t. Further reinforcing the validity of the assumption, the components of the first input vector, v1 , corresponding to the largest singular value, were essentially constant in time as compared to the components of the other input vectors, implying that v1 had converged to the time-invariant unit vector of period sensitivities as predicted. These results are demonstrated in Fig. 4, where σ1 dominates σ2 by an order of magnitude after about 7τ. Note that, in the SVD-based method determination of the sign of
S τ also requires direct perturbation of one of the parameters, which has already been performed for the results below. Since the SVD calculation derives S τ at a given point in time, it was possible to explore the time-dependence of the S τ computed using that method, as it was for the method of Larter (1983). Plots showing the time courses of the elements of S τ calculated using both methods are shown in Fig. 5 for parameter set A. These figures demonstrate that the period sensitivities calculated by the two methods have very different temporal behaviors. For the method of Larter (1983), the computed S τ was unbounded in time and median values had to be calculated (with the sole exception of vs in parameter set A). For the SVD-based method, the elements of S τ converged over time, and thus the value of the period sensitivity was taken to be that at large t. In the present study, the el-
Fig. 4. (a) First and second singular values for raw sensitivity matrix over time. Clearly, at large times, the first singular value dominates the second by several orders of magnitude. (b) Absolute values of the elements corresponding to ks of the first (v1 ) and second (v2 ) input vectors of the raw sensitivities over time. The element of the input vector v1 corresponding to the dominant singular value is essentially constant over time, while the element of the v2 vector shows complex variation.
670
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
Fig. 5. Period sensitivities for the parameters vs (a) and k2 (b) over time calculated using the two methods. Method 1 (Larter (1983)); Method 2 (SVD).
ements of S τ converged to their correct values after a time of approximately 24τ. In Fig. 6a, the values of the period sensitivities calculated using the two methods for parameter set A are plotted against each other. It is clear that the two approaches yield nearly the same values (similar results were obtained for parameter set B). To confirm the period sensitivity results for the two methods, period sensitivities were also computed using the finite difference method (Eq. (3)) for 1% perturbations in each of the parameters. Note that sensitivity behavior is generally a nonlinear characteristic, and the sensitivity analysis performed using the method of Larter (1983) or SVD are only linear analyses whose results will be valid locally in parameter space. For this reason, differences between the results of the sensitivity analysis and direct perturbation of the parameters may be expected. Nevertheless, very close agreement between period sensitivities calculated using the SVD-based method and those calculated using perturbations for parameter set A was observed (similar results were obtained for parameter set B), shown in Fig. 6b. To further validate the SVD method for computing period sensitivities, we also applied it to the Brusselator (Larter et al., 1984) and
Oregonator (Edelson & Thomas, 1981) model systems and obtained results that agreed with the literature (results not shown). It was also of interest to explore how the relative ordering of the period sensitivities varied between the two parameter sets. To facilitate comparison, we first computed scaled period sensitivities, S sτ , the components of which are given by: sτsj =
sτj pj τ
(35)
In Fig. 7a, the scaled sensitivities for the two parameter sets are plotted against each other. Note that the relative ordering of the parameters is largely invariant. This result is highly unexpected, given the nonlinear nature of the sensitivities and the local validity of the results. This result is suggestive of a conservation of sensitivities across the parameter sets. Finally, we wanted to investigate how the period of the oscillator would be altered by vectorial versus scalar perturbations. We first defined a scaled perturbation strength ps , which is a measure of the distance between the nominal parameter set p0 and the perturbed parameter set p:
Fig. 6. (a) Period sensitivities calculated using the SVD-based method are plotted against those calculated using the method of Larter (1983) for parameter set A. Note that the agreement is very good for both parameter sets. (b) Period sensitivities calculated from 1% parameter perturbations are plotted against those obtained using the SVD-based method for parameter set A. In general, the agreement is excellent.
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
671
Fig. 7. Scaled period sensitivities and oscillator robustness. (a) Scaled period sensitivities (circles) were calculated for the Goldbeter (1995) circadian rhythm oscillator, using parameter sets A and B. The line indicates perfect match. (b) Oscillator precision as a function of perturbation strength. Either single parameters showing the highest and second-highest absolute period sensitivities (squares and open circles, respectively) were perturbed, or all model parameters were modified in the direction of their period sensitivities (filled circles). Termination of the series indicates that stronger perturbations abolished oscillations altogether.
np pi − p0 2 i s p = p0i
(36)
i
The perturbed parameter vector was constructed as: pi = p0i (1 + ϕvpi )
(37)
where ϕ is a scaling factor, and vp is a np × 1 vector that defines the direction of the perturbations. For scalar perturbations, vp contained only single non-zero entries which were selected according to the relative period sensitivities of all parameters. When vectorial perturbations were investigated, vp was Sτ . For a given ps and vp , ϕ is given by: ps . ϕ=
vp T vp
(38)
Finally, oscillator precision was quantified by the scaled period deviation: τ(p) − τ(p0 ) . τ s = (39) τ(p ) 0
where τ s = 10−2 , for example, indicates a 1% change in period length relative to the unperturbed value. Results of the perturbation analysis for parameter set A are shown in Fig. 7b. Over a wide range of perturbation strengths, the perturbation in the direction of the period sensitivities clearly had a greater impact on the system period than a scalar perturbation to the most sensitive parameter of the same magnitude. This was contrasted, however, by the effect of high strength scalar perturbations that abolished oscillations completely. 4.1. Conclusions In the present work, we considered the sensitivity analysis of biological oscillators, focusing on the sensitivity of oscillator period to variations in model parameters. We discussed the current methods for calculating period sensitivities (Edelson
& Thomas, 1981; Larter, 1983) and their drawbacks in that the period sensitivities must be computed individually, care must be taken in selecting the states that are used and the times at which the period sensitivities are evaluated, and that the sensitivities may need to be integrated for at least one state (Larter, 1983). We described the development of a novel singular value decomposition (SVD)-based method for computing period sensitivities that offers several advantages over the current methods. The SVD-based method only requires that the period be sensitive to at least one of the parameters and simultaneously gives values for all the period sensitivities that converge asymptotically over time to their true values. For these reasons, the SVD-based method is easier to implement and interpret than the current methods. The drawbacks are that only the relative, not absolute, signs of the period sensitivities are obtained from the SVD-based method; and that it is generally more computationally intensive since it uses the sensitivities of all of the states. Like the previous methods (Edelson & Thomas, 1981; Larter, 1983), the SVD-based method also involves a transient during which the estimates of the period sensitivities converge to their true values. Determination of the absolute signs of the period sensitivities may be accomplished by at least one finite difference perturbation to the system. The convergence of the period sensitivities calculated using the SVD-based method is determined by how quickly the first term on the right hand side of Eq. (12) overwhelms the second term. Since the first term is dependent on the specific properties of the system being analyzed, there is no general rule for how large of t is required. For this reason, convergence of S τ is best assessed by inspection. We observed that, generally, this transient was longer for the SVD-based method than the previous methods (results not shown), although the results of the method of Larter (1983) were rapidly convoluted by the unbounded growth of the estimates. Nevertheless, we feel the SVD-based method is more straightforward than the current methods, since it is easier to check for convergence than it is to select appropriate times for evaluation of the
672
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
S τ or to address the unbounded growth of the computed result. As a case study, we performed a sensitivity analysis of a published model for circadian rhythms (Goldbeter, 1995). Our results showed an approximately constant relative ordering of period sensitivities for two parameter sets, suggesting that for this system, model structure, rather than location in parameter space, may be the key determinant of which parameters are most critical. Additionally, we performed a perturbation analysis of the system that demonstrated how the period sensitivity vector is the most sensitive or weak direction of the system. These results demonstrate how sensitivity analysis of computational models of biological oscillators may yield insights that are relevant to biological function. Given the importance of oscillatory systems in biology, and given the increased application of computational methods in biology, the importance of the analysis of computational models of biological oscillators can only be expected to grow. Sensitivity analysis is an important tool that may be used to gain understanding of mechanisms behind complex biological oscillators. The SVD-based method described in the present work simplifies the implementation and interpretation of sensitivity analysis of biological oscillators and thus may play a role in the further study of these important systems.
FJD acknowledges the support of an Alexander von Humboldt Research Fellowship and funding from the DARPA BioComp Initiative and the Institute for Collaborative Biotechnologies through grant DAAD19-03-D-0004 from the U.S. Army Research Office. DEZ acknowledges NIH training grant NIAAA5T32AA07463-15 and the University of Delaware Department of Chemical Engineering for funding. Appendix A The equations for the Goldbeter (1995) model are are as follows (Goldbeter, 1995): Kn dMP MP = vs n I n − vm dt Km + M P KI + P N dP0 P0 P1 = ks MP − v1 + v2 dt K1 + P 0 K2 + P 1
dP2 P1 P2 = v3 − v4 − k 1 P2 dt K3 + P 1 K4 + P 2 P2 + k2 PN − vd Kd + P 2 dPN = k1 P2 − k2 PN dt
Parameter
A
B
vs KI n vm Km ks v1 K1 v2 K2 v3 K3 v4 K4 vd Kd k1 k2
0.76 1 4 0.65 0.5 0.38 3.2 2 1.58 2 5 2 2.5 2 0.95 0.2 1.9 1.3
0.5 2 4 0.3 0.2 2 6 1.5 3 2 6 1.5 3 2 1.5 0.1 2 1
Values of the parameters for parameters sets A and B are given in Table 1.
References
Acknowledgments
dP1 P0 P1 P1 = v1 − v2 − v3 dt K1 + P 0 K2 + P 1 K3 + P 1 P2 + v4 K4 + P 2
Table 1 Parameter sets for the Goldbeter model, A (Goldbeter, 1995) and B (Gonze et al., 2002b)
(A.1)
Barkai, N., & Leibler, S. (2000). Circadian clocks limited by noise. Nature, 403, 267–268. Edelson, D., & Thomas, V. M. (1981). Sensitivity analysis of oscillating reactions 1: The period of the oregonator. J. Phys. Chem., 85 (11), 1555– 1558. Elowitz, M. B., Levine, A. J., Siggia, E. D., & Swain, P. S. (2002). Stochastic gene expression in a single cell. Science, 297, 1183–1186. Forger, D. B., & Peskin, C. S. (2003). A detailed predictive model of the mammalian circadian clock. Proc. Natl. Acad. Sci. USA, 100 (25), 14806–14811. Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22, 403–434. Goldbeter, A. (1995). A model for circadian oscillations in the Drosophila period protein PER. Proc. R. Soc. Lond. B., 261 (1362), 319–324. Goldbeter, A. (1996). Biochemical oscillations and cellular rhythms: The molecular bases of periodic and chaotic behaviour. (2nd ed.). Cambridge: Cambridge University Press. Goldbeter, A. (2002). Computational approaches to cellular rhythms. Nature, 420 (6912), 238–245. Gonze, D., Halloy, J., & Goldbeter, A. (2002). Deterministic versus stochastic models for circadian rhythms. J. Biol. Phys., 28, 637–653. Gonze, D., Halloy, J., & Goldbeter, A. (2002). Robustness of circadian rhythms with respect to molecular noise. Proc. Natl. Acad. Sci. USA, 99 (2), 673–678. Gonze, D., Halloy, J., Leloup, J. C., & Goldbeter, A. (2003). Stochastic models for circadian rhythms: effect of molecular noise on periodic and chaotic behaviour. C. R. Biol., 326 (2), 189–203. Hastings, M. H. (2000). Circadian clockworks: Two loops are better than one. Nat. Rev. Neurosci, 1, 143–146. Khalil, H. K. (1992). Nonlinear Systems. New York, N.Y.: Macmillan. Kholodenko, B. N., Demin, O. V., & Westerhoff, H. V. (1997). Control analysis of periodic phenomena in biological systems. J. Phys. Chem. B., 101 (11), 2070–2081.
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673 Larter, R. (1983). Sensitivity analysis of autonomous oscillators: Separation of secular terms and determination of structural stability. J. Phys. Chem., 87 (16), 3114–3121. Larter, R., Rabitz, H., & Kramer, M. (1984). Sensitivity analysis of limit cycles with application to the Brusselator. J. Chem. Phys., 80 (9), 4120– 4128. Leif, F., & Jorgensen, S. B. (2001). Estimation of kinetic parameters in a structured yeast model using regularisation. J. Biotechnol., 88 (3), 223– 237. Leloup, J. C., & Goldbeter, A. (1999). Chaos and birhythmicity in a model for circadian oscillations of the PER and TIM proteins in Drosophila. J. Theor. Biol., 198 (3), 445–459. Leloup, J. C., & Goldbeter, A. (2003). Toward a detailed computational model for the mammalian circadian clock. Proc. Natl. Acad. Sci. USA, 100 (12), 7051–7056. Levsky, J. M., Shenoy, S. M., Pezo, R. C., & Singer, R. H. (2002). Single-cell gene expression profiling. Science, 297 (5582), 836–840. Maly, T., & Petzold, L. R. (1996). Numerical methods and software for sensitivity analysis of differential–algebraic systems. Appl. Numer. Math., 20 (1/2), 57–79. Morohashi, M., Winn, A. E., Borisuk, M. T., Bolouri, H., Doyle, J., & Kitano, H. (2002). Robustness as a measure of plausibility in models of biochemical networks. J. Theor. Biol., 216 (1), 19–30. Novak, B., Pataki, Z., Ciliberto, A., & Tyson, J. J. (2001). Mathematical model of the cell division cycle of fission yeast. Chaos, 11 (1), 277–286. Rabitz, H., & Edelson, D. (1985). Numerical techniques for modelling and analysis of oscillating chemical reactions. In R. J. Field & M. Burger (Eds.), Oscillations and traveling waves in chemical systems (pp.193– 222). John Wiley and Sons. Rao, C. V., Wolf, D. M., & Arkin, A. P. (2002). Control, exploitation and tolerance of intracellular noise. Nature, 420 (6912), 231–237.
673
Savageau, M. A. (1971). Parameter sensitivity as a criterion for evaluating and comparing the performance of biochemical systems. Nature, 229 (5286), 542–544. Schlosser, P. M. (1994). Experimental design for parameter estimation through sensitivity analysis. J. Toxicol. Environ. Health, 43 (4), 495– 530. Shampine, L. F., & Reichelt, M. W. (1997). The MATLAB ODE suite. SIAM J. Sci. Comput., 18 (1), 1–22. Skogestad, S., & Postlethwaite, I. (2000). Multivariable Feedback Control. Chichester, UK: John Wiley and Sons. Stelling, J., & Gilles, E. D. (2001). Robustness vs. identifiability of regulatory modules?. Proc. 2nd Intl. Conf. Systems Biology 181–190. Sveiczer, A., Tyson, J. J., & Novak, B. (2001). A stochastic, molecular model of the fission yeast cell cycle: role of the nucleocytoplasmic ratio in cycle time regulation. Biophys. Chem., 92 (1/2), 1–15. Tomovi´c, R., & Vukobratovi´c, M. (1972). General Sensitivity Theory. New York, N.Y.: American Elsevier. Tyson, J. J., Csikasz-Nagy, A., & Novak, B. (2002). The dynamics of cell cycle regulation. Bioessays, 24 (12), 1095–1109. Vilar, J. M., Kueh, H. Y., Barkai, N., & Leibler, S. (2002). Mechanisms of noise–resistance in genetic oscillators. Proc. Natl. Acad. Sci. USA, 99 (30), 5988–5992. Wang, S. H., & Lee, T. F. (1996). Interpretation of non–minimum phase system with extension to time–varying case. Electron. Lett., 32 (3), 270– 272. Wang, S. H., & Zachery, R. (1996). Singular value decomposition of system input–output matrix and its symmetry property. Comput. Electron. Eng., 22 (3), 231–234. Zak, D. E., Doyle, F. J., Vlachos, D. G., & Schwaber, J. S. (2001). Stochastic kinetic analysis of transcriptional feedback models for circadian rhythms. Proc. 40th IEEE Conf. Decision & Control 849–854.
Sensitivity analysis of oscillatory (bio)chemical systems Daniel E. Zaka , J¨org Stellingb , Francis J. Doyle IIIc,∗ b
a Department of Chemical Engineering, University of Delaware, Newark, DE 19716, USA Max Planck Institute for Dynamics of Complex Technical Systems, D-39106 Magdeburg, Germany c Department of Chemical Engineering, University of California, Santa Barbara, CA 93106, USA
Received 22 September 2003; received in revised form 8 January 2004 Available online 7 October 2004
Abstract From the cell cycle to circadian rhythms, oscillatory processes are fundamental to biology. Emerging from nonlinear dynamical interactions, oscillatory mechanisms are best understood through mathematical modeling. Ordinary differential equations (ODEs) are one framework in which the complex interactions giving rise to biological oscillations may be modeled. Key to ODE models are the model parameters that determine whether or not oscillations will occur, and the period and amplitude of the oscillations when they do. Sensitivity analysis is a means to acquire insight about the importance of the model parameters. Sensitivity analysis of oscillatory systems provides unique challenges and must be addressed carefully. In the present study, we describe a method for determining the sensitivity of the period to the model parameters that is straightforward to implement and interpret. We apply this method to a model for circadian rhythms, and obtain results suggesting a link between network structure and parameter sensitivity. © 2004 Elsevier Ltd. All rights reserved. Keywords: Circadian rhythm; Mathematical modeling; Sensitivity analysis
1. Introduction Rhythmic processes are fundamental to biology, as is revealed by literature searches on two key oscillatory biological processes, ‘circadian rhythm’ and ‘cell cycle’, that yielded 43,057 and 61,556 articles, respectively (http://www.ncbi. nlm.nih.gov: December 22, 2003). Circadian rhythms are variations in activity (physiological and cellular) with nearly twenty-four hour periods that, for the fly and the mouse, arise from cellular genetic networks containing delayed feedback mechanisms (Hastings, 2000). The cell cycle is the process by which a single cell divides to become two, and is comprised of a complex dynamic nonlinear network of protein interactions (Tyson, Csikasz-Nagy, & Novak, 2002). The oscillations in these and other cellular processes arise from biochemical reactions, and thus general principles that have been developed for chemical oscillators apply to biological oscillators. One ∗
Corresponding author. E-mail address: [email protected] (F.J. Doyle).
0098-1354/$ – see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2004.08.021
principle is that computational approaches are necessary to fully understand and explain their behavior (Rabitz & Edelson, 1985; Goldbeter, 1996; Goldbeter, 2002). The approaches employed for the mathematical analysis of biological oscillators have paralleled those used for chemical oscillators, where both ordinary differential equation (ODE) formulations, derived from mass action kinetics, and discrete stochastic formulations (Gillespie, 1976) have been employed, with the ODE approaches being predominant (Rabitz & Edelson, 1985). Interest in stochastic models of biological oscillators (Barkai & Leibler, 2000; Gonze, Halloy, & Goldbeter, 2002b; Vilar, Kueh, Barkai, & Leibler, 2002; Zak, Doyle, Vlachos, & Schwaber, 2001), and of biological systems in general, is growing, however, and has been strengthened by the publication of several elegant experimental studies that demonstrated the stochastic nature of some biochemical reactions at the single-cell level (Levsky, Shenoy, Pezo, & Singer, 2002; Elowitz, Levine, Siggia, & Swain, 2002), reviewed in (Rao, Wolf, & Arkin, 2002). In spite of this, new models of biological oscillators in the ODE framework continue to be developed (Leloup & Goldbeter,
664
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
Nomenclature A B f p p0 S(t) Sτ S sτ S c (t) sij sτj sτsj v vp x
ns × ns state Jacobian matrix ns × np parameter Jacobian matrix vector of state derivatives of size ns vector of parameters of size np nominal parameter values ns × np sensitivity matrix period sensitivity vector scaled period sensitivity vector ns × np cleaned out sensitivity matrix sensitivity of the ith state to the jth parameter sensitivity of the period to the jth parameter scaled sensitivity of the period to the jth parameter input direction from SVD perturbation direction vector state vector of size ns
Greek letters ps scaled perturbation strength τ s scaled period deviation (oscillator precision) ϕ perturbation scaling factor σ singular value τ period of the oscillator
mechanistic insight and to design informative experiments. Sensitivity analysis is one technique to investigate the importance of parameters. Sensitivity analysis has been applied in a few cases to the analysis of biological systems, with the objectives of mechanism discrimination on the basis of model sensitivity (Morohashi et al., 2002; Savageau, 1971), experimental design and parameter estimation (Leif & Jorgensen, 2001; Schlosser, 1994), and the relationship between sensitivity and identifiability (Stelling & Gilles, 2001). In the present work, we consider the parametric sensitivity analysis of ODE models of biological oscillators, where the parameters determine both the existence of and characteristics of the oscillations. Since the period of biological oscillators is a key aspect of their physiological significance (for example, the time-keeping nature of circadian rhythms), we specifically consider the parametric sensitivity of the period. Direct application of sensitivity analysis to oscillating systems gives rise to secular terms (Larter, 1983), and thus must be done with care. In the present work, we first discuss the standard methods for the sensitivity analysis of oscillatory systems, and then present a novel method based on singular value decomposition (SVD) that is easier to interpret and implement than the common methods. We conclude with a case study of an ODE model for circadian rhythms. 2. Sensitivity analysis of oscillatory systems
2003; Forger & Peskin, 2003), and stochastic and deterministic simulations are often used as complementary, rather than competing, methods by individual groups (Gonze, Halloy, & Goldbeter, 2002a; Gonze, Halloy, Leloup, & Goldbeter, 2003; Leloup & Goldbeter, 1999; Novak, Pataki, Ciliberto, & Tyson, 2001; Sveiczer, Tyson, & Novak, 2001; Vilar et al., 2002). Generally, the ODE framework is used because of its ease of simulation and analysis with available software tools (Maly & Petzold, 1996; Shampine & Reichelt, 1997), while the stochastic framework is used to explore fluctuations. Interestingly, there have been a number of studies of biological oscillators where predictions made using the ODE framework have held up in a stochastic framework, with the effect of stochastic noise largely being the widening of the limit cycle (Gonze, Halloy, & Goldbeter, 2002a; Gonze et al., 2003). Thus, while experimental results have demonstrated that the stochastic nature of some biological reactions cannot be neglected, ODEs, and tools to analyze them, are likely to continue to be of importance to the computational biology community. Given that oscillating chemical systems typically consist of larger numbers of molecules and potentially faster reactions, we also expect ODEs to continue to play an important role in the study of chemical oscillators. ODE models of biological processes often involve many parameters, and the importance of these parameters in determining system behavior must be assessed in order to gain
In the present section, we first describe basic principles of sensitivity analysis that are applicable to any ODE system. This is followed by a presentation of aspects of sensitivity analysis that are specific to oscillatory systems. 2.1. General principles of sensitivity analysis Consider the ODE system: x˙ = f (x, p)
(1)
where x(t) is a vector of ns states, p is a vector of np model parameters, and f is a column vector of the state time derivatives. Assuming that a solution to Eq. (1) exists, the sensitivity matrix of the system, S(t), that describes how variations in the parameters near the local point in parameter space, p0 , influence the state trajectories, may be defined: ∂x (2) S(t) ≡ ∂p x=x(t,p0 ),p=p0 where S(t) is composed of individual sensitivities of each state to each parameter (sij ). The simplest way to calculate S(t) is by finite differences. For a single parameter pj , x(t, pj + pj ) − x(t, pj ) ∂x ≈ (3) ∂pj x=x(t,p ) pj 0
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
While Eq. (3) may be used to estimate the sensitivities, it is computationally tedious and frequently inaccurate, given that it may lead to numerical instabilities (Rabitz & Edelson, 1985). An alternative approach is to differentiate Eq. (1) with respect to the parameter vector p, giving (Khalil, 1992): S˙ = A(t, p0 )S + B(t, p0 )
2.2. Oscillatory systems The present work is concerned with systems that are periodic in time, for which: x(t + τ) = x(t)
(6)
where τ is the period of the oscillation. The objective is to determine how τ depends on the model parameters, p, specifically to determine the np × 1 vector of period sensitivities, Sτ : T ∂τ ∂τ ∂τ Sτ ≡ (7) , ,···, ∂p1 ∂p2 ∂pnp S τ is composed of individual parameter sensitivities, sτj , and is a constant vector. Given that x(t) is periodic, it is possible to express each of the states x1 , x2 , . . . xns as Fourier series (Larter, 1983):
n=0
ani cos
2nπt 2nπt + bni sin τ τ
(8)
where ani and bni are Fourier coefficients and are functions of parameters. The time derivative x˙ i can also be calculated: ∞ 2π 2nπt 2nπt −nani sin + nbni cos τ τ τ n=0
n=0
∂ani 2nπt ∂bi 2nπt cos + n sin ∂pj τ ∂pj τ
(10)
Using Eqs. 8 and 9, Eq. (10) can be rearranged to: (5)
The sensitivity matrix S(t) as defined by Eq. (4) is a linear time-varying (LTV) system. By symbolically calculating the state and parameter Jacobians (A and B, respectively) it is possible to integrate x(t) and S(t) simultaneously, given the nominal parameter values, p0 , and the initial conditions for the states, x0 . Unless the initial conditions of the states depend on the parameters, the initial conditions for the sensitivities are zero. This method is known as the direct method (Rabitz & Edelson, 1985), and involves the integration of a system of dimension ns × (np + 1). When np > ns , Eq. (5) may be more efficiently integrated through the use of a Green’s function technique (Rabitz & Edelson, 1985).
∞
∞ n=0
∂f ∂x x=x(t,p0 ),p=p0 ∂f B(t, p0 ) = ∂p x=x(t,p0 ),p=p0
x˙ i (t) =
∞ 2πt ∂τ 2nπt 2nπt i i sij = 2 nan sin − nbn cos τ τ τ ∂pj
(4)
A(t, p0 ) =
xi (t) =
Finally, the sensitivity sij may be calculated, keeping in mind that τ depends on p:
+
where
665
(9)
t sij = − sτj fi + τ
∂xi ∂pj
(11) τ
or in vector notation: t S = − f S Tτ + S c τ
(12)
where S c (t) is the ns × np cleaned out sensitivity matrix (Tomovi´c & Vukobratovi´c, 1972), composed of (∂xi /∂pj )τ , and is periodic in time, as can be seen in Eq. (10). S c (t) contains information about how variations in the parameters affect the shape of the trajectory (Larter, 1983). The above Fourier series-based derivation was given by Larter (1983), but Eq. (11) was originally derived by Tomovi´c and Vukobratovi´c (1972). Eq. (12) clearly demonstrates that, because of the secular term, the raw sensitivity coefficients sij are unbounded in time for parameters with sτj = 0. A similar observation has been made in the metabolic control analysis of autonomous oscillators (Kholodenko, Demin, & Westerhoff, 1997). The secular term arises because the computation of the sensitivities involves a nonuniformly valid expansion of a periodic system (Larter, 1983). Through the Linstedt–Poincar´e method, Larter (1983) has shown that the periodic cleaned out sensitivity matrix, S c (t) is actually the first correction term in a uniformly valid expansion of the system around the nominal parameter set (Larter, 1983). To determine either how variation in the parameters influences the oscillator period (S τ ) or the shape of its trajectory (S c (t)), the raw sensitivity matrix (S) must be broken into its component terms. This is generally accomplished by calculating S τ by another means, and then applying a correction to S(t) (Larter, 1983). Two methods for computing S τ have been described in the literature, and these will be considered in the following section. 2.3. Current methods for determining period sensitivities Edelson and Thomas (1981) have described one approach for calculating S τ . Differentiating Eq. (6) for a particular state, xi , with respect to a particular parameter,
666
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
pj , gives: dxi (t, p) dxi (t + τ, p) = dpj dpj ∂xi ∂t ∂xi ∂xi ∂(t + τ) ∂xi + = + ∂t ∂pj ∂pj t ∂(t + τ) ∂pj ∂pj t+τ ∂xi ∂xi ∂τ ∂xi = + ∂pj t ∂t ∂pj t ∂pj t+τ
Eq. (16) may be simplified by combining with Eq. (11). Firstly: sτj ∂xi (17) txi − xi dt + sij dt = − τ0 ∂pj τ (13)
Rearranging gives:
=
Sij (t) − Sij (t + τ) fi (t)
(14)
In practice, Eq. (14) additionally will contain a transient term that decays over time to yield the correct sτj , even if the initial conditions lie on the limit cycle (Rabitz & Edelson, 1985). This transient results from the fact that these initial conditions are actually dependent on the model parameters (Rabitz & Edelson, 1985). It must be pointed out that computation of sτj by Eq. (14) involves division by the time derivative of xi , fi , that must be zero at some points along the trajectory because the system is oscillatory. For this reason, period sensitivities computed using this method are not defined at all points along the cycle. As described by Edelson and Thomas (1981) themselves, care also must be taken in choosing which state is used to compute the period sensitivities, to avoid both division by small time derivatives and inaccuracies in the difference in the numerator. For the Oregonator system they considered, two out of the three states of the dynamic system proved to be undesirable for computing period sensitivities (Edelson & Thomas, 1981). An alternative integral approach to computing period sensitivities that removes the need to divide by time derivatives (but not does not resolve the issue of initial transients) has been given by Larter (1983): t1 +τ sτj = −
t1
sij dt −
t2 +τ t2
sij dt
(15)
xi (t1 ) − xi (t2 )
Eq. (15) will not exhibit instabilities as long as t1 and t2 are chosen so that xi (t1 ) − xi (t2 ) is not small. This improvement in numerical performance is traded off against an increased demand on the accuracy of the period (τ) used in the integrations, however. This aspect of Eq. (15) arises through the numerator, and is explained in the following derivation. Suppose that the period, τ, used for the integrations in Eq. (15) is actually the true period of the oscillator, τ0 , plus some small error . This gives:
tk +τ
tk
sij dt =
tk +τ0
tk
sij dt +
tk tk +τ0 tk
(∂xi /∂pj )t − (∂xi /∂pj )t+τ = (∂xi /∂t)t
sτj
Since xi and (∂xi /∂pj )τ are both periodic, the following are true for any time tk : tk +τ0 xi dt = Cx
tk +τ0 + tk +τ0
sij dt
(16)
∂xi ∂pj
(18) τ
dt = Cτ
where both Cx and Cτ are constants. The following are reasonable approximations when is small: tk +τ0 + xi dt ≈ Cx tk tk +τ0 + ∂x i
tk
∂pj
τ
dt ≈ Cτ
(19)
xi (tk ≤ t ≤ tk + ) ≈ fi (tk )t + b(tk ) where fi (tk ) and b(tk ) are from the local linearization of xi in time and are periodic. Combining Eq. (16) with Eqs. 17–19, and neglecting all terms of order and higher, gives: tk +τ sτj sτj sij dt ≈ −sτj xi (tk )+ Cx + Cτ − fi (tk )tk (20) τ τ0 0 tk Substituting Eq. (20) into Eq. (15) gives: sτcalc = sτj j
(xi (t1 ) − xi (t2 )) + (/τ0 )(fi (t1 )t1 − fi (t2 )t2 ) (21) xi (t1 ) − xi (t2 )
In general, fi (t1 ) and fi (t2 ) will differ, and thus the error will grow unbounded with increasing t1 and t2 . The rate at which the error grows can be slowed by choosing t1 and t2 so that fi (t1 ) ≈ fi (t2 ). To summarize, the current methods for computing S τ involve calculating sτj for each parameter individually, and require careful selection of the states and times at which the calculations are performed. For the method of Edelson & Thomas (1981), t must be chosen so that fi is not close to zero. For the method of Larter (1983), the integrated sensitivities for at least one state xi (t) must be computed, and t1 and t2 must be chosen to maximize the difference between xi (t1 ) and xi (t2 ), while minimizing the difference between fi (t1 ) and fi (t2 ). A method that can determine all of the period sensitivities at once and without concerns of numerical stability is desirable and is the focus of the following section. 3. Singular value decomposition approach to period sensitivity The novel method for determining period sensitivities that is proposed in the present work is based on the singular value decomposition (SVD) of the raw sensitivity matrix S(t).
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
667
SVD has seen widespread application in multivariable feedback control, to identify strong and weak input directions, as well as to quantify interactions between inputs and outputs (Skogestad & Postlethwaite, 2000). SVD has also been applied to linear and LTV discrete time systems to gain a better understanding of non-minimum phase systems in the time domain (Wang & Lee, 1996) and the reachable subspace of system outputs (Wang & Zachery, 1996). The following derivation is based on the observation that, at large times, the first term on the right-hand side of Eq. (12) will dominate the periodic term S c , giving:
Postlethwaite, 2000), the single non-zero eigenvalue of S τ S Tτ is simply S Tτ S τ . This gives:
t S = − f S Tτ τ
v = αS τ
(22)
The value of t that constitutes sufficiently large will depend on the characteristics of the system being analyzed, specifically the magnitudes of τ, f , and S τ . To be valid, Eq. (22) also requires f = 0, and S τ = 0. Given that x(t) is periodic, the first condition obviously cannot be guaranteed for all points along the trajectory for all components of f . However, the periodicity of the system also guarantees that all components of f will never be zero simultaneously, because that would imply steady-state behavior. Except for some pathological systems (Larter, Rabitz, & Kramer, 1984), τ will be sensitive to at least one parameter, and it is a valid assumption that S τ = 0. The rare cases where it is not will be obvious from inspection of S(t), which will be purely periodic. The SVD of S(t) gives: S = UΣV T
(23)
where Σ is an ns × np diagonal matrix of non-negative singular values, σk , that are the square roots of the eigenvalues of S T S or SS T . V is an np × np matrix whose columns are unit eigenvectors of S T S. U is an ns × ns matrix whose columns are unit eigenvectors of SS T . The vectors which comprise U and V are termed the input and output vectors of S given that: Svk = σk uk
(24)
Thus, given an input along direction vk , the output will be in the direction uk with gain σk (Skogestad & Postlethwaite, 2000). Given that SVD of S(t) requires operations on the matrix S T S, it is convenient to simplify it using Eq. (22): 2 t2 T T 2t (25) S f f S = φ S τ S Tτ τ τ τ2 τ2
where φ(t) = f T f = 0 (since φ ≡ 0 would imply steadystate). It follows that:
STS =
eig(S T S) = φ2
t2 eig(S τ S Tτ ) τ2
(26)
where the expression eig indicates eigenvalues. Since, for any matrix X, the non-zero eigenvalues of XT X are equivalent to the eigenvalues of XXT (Skogestad &
t2 T (27) S Sτ τ2 τ To find the eigenvector v that corresponds to λ, and therefore the most sensitive input direction, the following condition must be satisfied:
eig(S T S) = λ = φ2
S T Sv = λv
(28)
The solution is: (29)
where α is an arbitrary scalar. The solution is readily verified: S T SαS τ = αφ2
t2 S τ S Tτ S τ = αλS τ = λv τ2
(30)
Thus, singular value decomposition of the raw sensitivity matrix at large t will yield a scalar multiple of the vector of period sensitivities, S τ , as the input vector corresponding to the only nonzero singular value. The determination of the magnitudes of the period sensitivities from the SVD now follows. Since eigenvectors are usually normalized to be unit vectors, the scalar multiple α will be:
−1 T α=± Sτ Sτ (31) Thus the input vector corresponding to the largest singular value will be: v = ±
Sτ
(32)
S Tτ S τ
where the relative, but not absolute, signs of the period sensitivities are known. The magnitudes of the period sensitivities can be found using the singular value. Given that:
√ t σ= λ=φ S Tτ S τ (33) τ it follows that: στ Sτ = ± v φt
(34)
In conclusion, the magnitudes and relative signs of the period sensitivities of an oscillatory system may be determined using the SVD of the raw sensitivity matrix at large t, when f = 0, and when τ is sensitive to at least one parameter. This result has an interesting interpretation in terms of the classical application of SVD in that the vector of period sensitivities corresponds to the most sensitive input direction of the raw sensitivity matrix. What constitutes sufficiently large for t will depend on the specific system analyzed, and thus should be determined by considering the convergence of the calculated values of S τ . Also, when Eq. (22) is strictly true, there will only be one non-zero singular value. Eq. (22) is an approximation to a system truly described by Eq. (12), however, and
668
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
Fig. 1. Raw sensitivities of the Goldbeter (1995) model, parameter set A, displaying unbounded growth linear in time as predicted by Eq. (12).
thus there will be several non-zero singular values. At large t, the first singular value, σ1 will dominate the others by several orders of magnitude, and under these conditions σ1 and v1 may be substituted for σ and v in the above expressions. The absolute sign of S τ also must be calibrated, which may be accomplished by a small perturbation to one of the parameters and determining if it increases or decreases the period. Overall, the present method has the advantages over the previous methods in that less care and human intervention is required in selecting states and times for the computations, and that the period sensitivity estimates are guaranteed to converge to their true values. These advantages come at the cost of increased computational requirements, however. 4. Case study: period sensitivity of Goldbeter (1995) circadian rhythm model The present section is concerned with the computation of period sensitivities from the raw sensitivities for the Goldbeter (1995) model for intracellular circadian rhythms.
As mentioned above, circadian rhythms are periodic daily variations of activity that arise from intracellular feedback mechanisms (Hastings, 2000). The Goldbeter (1995) model is a 5-state system (ns = 5) with 18 parameters (np = 18). Two sets of parameter values have been given, with set A being from the original paper (Goldbeter, 1995), yielding a period of 23.7 h. Parameter set B is from a recent publication (Gonze et al., 2002b) and yields a period of 22.2 h. The model equations and parameters are provided in the appendix. Oscillations arise in the system from the general principle of delayed negative feedback. The raw sensitivities for the system were calculated using Eq. (4) by simultaneous integration of the original model and the sensitivity equations (the direct method) in MATLAB (The Mathworks, Inc.) using the function ode15s (Shampine & Reichelt, 1997). The sensitivities were unbounded in time, as predicted by Eq. (12), shown for parameter set A in Fig. 1. The method of Larter (1983) (Eq. (15)) was first used to compute S τ for the Goldbeter (1995) model. In the present calculations, t1 was varied continuously and t2 was chosen to maximize the difference |xi (t1 ) − xi (t2 )| within one τ of time, (t2 < t1 + τ), ensuring that the denominator in Eq. (15) did not approach zero. Although Eq. (15) predicts that the period sensitivities should be constant over time (after the decay of initial transients), they were not, as shown in Fig. 2 for parameter set A. For the parameter vs , for example (Fig. 2a), the period sensitivity showed damped oscillations, eventually settling down to a constant value. This parameter was the exception, as for the remainder of the parameters for parameter set A, and for all of the parameters for parameter set B, the computed period sensitivities grew unbounded in time. This is shown for parameter k2 in Fig. 2b. This unbounded growth was due to the amplification of small errors in the value of τ that was used in the integrations, as discussed above. To obtain actual values of S τ , median values were computed from the time courses. Given values for S τ and Eq. (12), it was possible to calculate the cleaned out sensitivities S c (t). A representative plot
Fig. 2. Period sensitivities to (a) vs (∂τ/∂vs ) and (b) k2 (∂τ/∂k2 ) calculated using the method of Larter (1983).
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
669
Fig. 3. Raw and cleaned out sensitivities (∂MP /∂ks ) for 5-state Goldbeter model, parameter set A. Note that the cleaned out sensitivity is periodic in time as expected.
of a periodic cleaned out sensitivity with its corresponding unbounded raw sensitivity is shown in Fig. 3. The SVD approach also was used to compute the S τ for the Goldbeter model. As predicted by the derivations, at large t, the first singular value, σ1 , dominated the others by several orders of magnitude, demonstrating that the assumption made in deriving Eq. (22) was valid at large t. Further reinforcing the validity of the assumption, the components of the first input vector, v1 , corresponding to the largest singular value, were essentially constant in time as compared to the components of the other input vectors, implying that v1 had converged to the time-invariant unit vector of period sensitivities as predicted. These results are demonstrated in Fig. 4, where σ1 dominates σ2 by an order of magnitude after about 7τ. Note that, in the SVD-based method determination of the sign of
S τ also requires direct perturbation of one of the parameters, which has already been performed for the results below. Since the SVD calculation derives S τ at a given point in time, it was possible to explore the time-dependence of the S τ computed using that method, as it was for the method of Larter (1983). Plots showing the time courses of the elements of S τ calculated using both methods are shown in Fig. 5 for parameter set A. These figures demonstrate that the period sensitivities calculated by the two methods have very different temporal behaviors. For the method of Larter (1983), the computed S τ was unbounded in time and median values had to be calculated (with the sole exception of vs in parameter set A). For the SVD-based method, the elements of S τ converged over time, and thus the value of the period sensitivity was taken to be that at large t. In the present study, the el-
Fig. 4. (a) First and second singular values for raw sensitivity matrix over time. Clearly, at large times, the first singular value dominates the second by several orders of magnitude. (b) Absolute values of the elements corresponding to ks of the first (v1 ) and second (v2 ) input vectors of the raw sensitivities over time. The element of the input vector v1 corresponding to the dominant singular value is essentially constant over time, while the element of the v2 vector shows complex variation.
670
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
Fig. 5. Period sensitivities for the parameters vs (a) and k2 (b) over time calculated using the two methods. Method 1 (Larter (1983)); Method 2 (SVD).
ements of S τ converged to their correct values after a time of approximately 24τ. In Fig. 6a, the values of the period sensitivities calculated using the two methods for parameter set A are plotted against each other. It is clear that the two approaches yield nearly the same values (similar results were obtained for parameter set B). To confirm the period sensitivity results for the two methods, period sensitivities were also computed using the finite difference method (Eq. (3)) for 1% perturbations in each of the parameters. Note that sensitivity behavior is generally a nonlinear characteristic, and the sensitivity analysis performed using the method of Larter (1983) or SVD are only linear analyses whose results will be valid locally in parameter space. For this reason, differences between the results of the sensitivity analysis and direct perturbation of the parameters may be expected. Nevertheless, very close agreement between period sensitivities calculated using the SVD-based method and those calculated using perturbations for parameter set A was observed (similar results were obtained for parameter set B), shown in Fig. 6b. To further validate the SVD method for computing period sensitivities, we also applied it to the Brusselator (Larter et al., 1984) and
Oregonator (Edelson & Thomas, 1981) model systems and obtained results that agreed with the literature (results not shown). It was also of interest to explore how the relative ordering of the period sensitivities varied between the two parameter sets. To facilitate comparison, we first computed scaled period sensitivities, S sτ , the components of which are given by: sτsj =
sτj pj τ
(35)
In Fig. 7a, the scaled sensitivities for the two parameter sets are plotted against each other. Note that the relative ordering of the parameters is largely invariant. This result is highly unexpected, given the nonlinear nature of the sensitivities and the local validity of the results. This result is suggestive of a conservation of sensitivities across the parameter sets. Finally, we wanted to investigate how the period of the oscillator would be altered by vectorial versus scalar perturbations. We first defined a scaled perturbation strength ps , which is a measure of the distance between the nominal parameter set p0 and the perturbed parameter set p:
Fig. 6. (a) Period sensitivities calculated using the SVD-based method are plotted against those calculated using the method of Larter (1983) for parameter set A. Note that the agreement is very good for both parameter sets. (b) Period sensitivities calculated from 1% parameter perturbations are plotted against those obtained using the SVD-based method for parameter set A. In general, the agreement is excellent.
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
671
Fig. 7. Scaled period sensitivities and oscillator robustness. (a) Scaled period sensitivities (circles) were calculated for the Goldbeter (1995) circadian rhythm oscillator, using parameter sets A and B. The line indicates perfect match. (b) Oscillator precision as a function of perturbation strength. Either single parameters showing the highest and second-highest absolute period sensitivities (squares and open circles, respectively) were perturbed, or all model parameters were modified in the direction of their period sensitivities (filled circles). Termination of the series indicates that stronger perturbations abolished oscillations altogether.
np pi − p0 2 i s p = p0i
(36)
i
The perturbed parameter vector was constructed as: pi = p0i (1 + ϕvpi )
(37)
where ϕ is a scaling factor, and vp is a np × 1 vector that defines the direction of the perturbations. For scalar perturbations, vp contained only single non-zero entries which were selected according to the relative period sensitivities of all parameters. When vectorial perturbations were investigated, vp was Sτ . For a given ps and vp , ϕ is given by: ps . ϕ=
vp T vp
(38)
Finally, oscillator precision was quantified by the scaled period deviation: τ(p) − τ(p0 ) . τ s = (39) τ(p ) 0
where τ s = 10−2 , for example, indicates a 1% change in period length relative to the unperturbed value. Results of the perturbation analysis for parameter set A are shown in Fig. 7b. Over a wide range of perturbation strengths, the perturbation in the direction of the period sensitivities clearly had a greater impact on the system period than a scalar perturbation to the most sensitive parameter of the same magnitude. This was contrasted, however, by the effect of high strength scalar perturbations that abolished oscillations completely. 4.1. Conclusions In the present work, we considered the sensitivity analysis of biological oscillators, focusing on the sensitivity of oscillator period to variations in model parameters. We discussed the current methods for calculating period sensitivities (Edelson
& Thomas, 1981; Larter, 1983) and their drawbacks in that the period sensitivities must be computed individually, care must be taken in selecting the states that are used and the times at which the period sensitivities are evaluated, and that the sensitivities may need to be integrated for at least one state (Larter, 1983). We described the development of a novel singular value decomposition (SVD)-based method for computing period sensitivities that offers several advantages over the current methods. The SVD-based method only requires that the period be sensitive to at least one of the parameters and simultaneously gives values for all the period sensitivities that converge asymptotically over time to their true values. For these reasons, the SVD-based method is easier to implement and interpret than the current methods. The drawbacks are that only the relative, not absolute, signs of the period sensitivities are obtained from the SVD-based method; and that it is generally more computationally intensive since it uses the sensitivities of all of the states. Like the previous methods (Edelson & Thomas, 1981; Larter, 1983), the SVD-based method also involves a transient during which the estimates of the period sensitivities converge to their true values. Determination of the absolute signs of the period sensitivities may be accomplished by at least one finite difference perturbation to the system. The convergence of the period sensitivities calculated using the SVD-based method is determined by how quickly the first term on the right hand side of Eq. (12) overwhelms the second term. Since the first term is dependent on the specific properties of the system being analyzed, there is no general rule for how large of t is required. For this reason, convergence of S τ is best assessed by inspection. We observed that, generally, this transient was longer for the SVD-based method than the previous methods (results not shown), although the results of the method of Larter (1983) were rapidly convoluted by the unbounded growth of the estimates. Nevertheless, we feel the SVD-based method is more straightforward than the current methods, since it is easier to check for convergence than it is to select appropriate times for evaluation of the
672
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673
S τ or to address the unbounded growth of the computed result. As a case study, we performed a sensitivity analysis of a published model for circadian rhythms (Goldbeter, 1995). Our results showed an approximately constant relative ordering of period sensitivities for two parameter sets, suggesting that for this system, model structure, rather than location in parameter space, may be the key determinant of which parameters are most critical. Additionally, we performed a perturbation analysis of the system that demonstrated how the period sensitivity vector is the most sensitive or weak direction of the system. These results demonstrate how sensitivity analysis of computational models of biological oscillators may yield insights that are relevant to biological function. Given the importance of oscillatory systems in biology, and given the increased application of computational methods in biology, the importance of the analysis of computational models of biological oscillators can only be expected to grow. Sensitivity analysis is an important tool that may be used to gain understanding of mechanisms behind complex biological oscillators. The SVD-based method described in the present work simplifies the implementation and interpretation of sensitivity analysis of biological oscillators and thus may play a role in the further study of these important systems.
FJD acknowledges the support of an Alexander von Humboldt Research Fellowship and funding from the DARPA BioComp Initiative and the Institute for Collaborative Biotechnologies through grant DAAD19-03-D-0004 from the U.S. Army Research Office. DEZ acknowledges NIH training grant NIAAA5T32AA07463-15 and the University of Delaware Department of Chemical Engineering for funding. Appendix A The equations for the Goldbeter (1995) model are are as follows (Goldbeter, 1995): Kn dMP MP = vs n I n − vm dt Km + M P KI + P N dP0 P0 P1 = ks MP − v1 + v2 dt K1 + P 0 K2 + P 1
dP2 P1 P2 = v3 − v4 − k 1 P2 dt K3 + P 1 K4 + P 2 P2 + k2 PN − vd Kd + P 2 dPN = k1 P2 − k2 PN dt
Parameter
A
B
vs KI n vm Km ks v1 K1 v2 K2 v3 K3 v4 K4 vd Kd k1 k2
0.76 1 4 0.65 0.5 0.38 3.2 2 1.58 2 5 2 2.5 2 0.95 0.2 1.9 1.3
0.5 2 4 0.3 0.2 2 6 1.5 3 2 6 1.5 3 2 1.5 0.1 2 1
Values of the parameters for parameters sets A and B are given in Table 1.
References
Acknowledgments
dP1 P0 P1 P1 = v1 − v2 − v3 dt K1 + P 0 K2 + P 1 K3 + P 1 P2 + v4 K4 + P 2
Table 1 Parameter sets for the Goldbeter model, A (Goldbeter, 1995) and B (Gonze et al., 2002b)
(A.1)
Barkai, N., & Leibler, S. (2000). Circadian clocks limited by noise. Nature, 403, 267–268. Edelson, D., & Thomas, V. M. (1981). Sensitivity analysis of oscillating reactions 1: The period of the oregonator. J. Phys. Chem., 85 (11), 1555– 1558. Elowitz, M. B., Levine, A. J., Siggia, E. D., & Swain, P. S. (2002). Stochastic gene expression in a single cell. Science, 297, 1183–1186. Forger, D. B., & Peskin, C. S. (2003). A detailed predictive model of the mammalian circadian clock. Proc. Natl. Acad. Sci. USA, 100 (25), 14806–14811. Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22, 403–434. Goldbeter, A. (1995). A model for circadian oscillations in the Drosophila period protein PER. Proc. R. Soc. Lond. B., 261 (1362), 319–324. Goldbeter, A. (1996). Biochemical oscillations and cellular rhythms: The molecular bases of periodic and chaotic behaviour. (2nd ed.). Cambridge: Cambridge University Press. Goldbeter, A. (2002). Computational approaches to cellular rhythms. Nature, 420 (6912), 238–245. Gonze, D., Halloy, J., & Goldbeter, A. (2002). Deterministic versus stochastic models for circadian rhythms. J. Biol. Phys., 28, 637–653. Gonze, D., Halloy, J., & Goldbeter, A. (2002). Robustness of circadian rhythms with respect to molecular noise. Proc. Natl. Acad. Sci. USA, 99 (2), 673–678. Gonze, D., Halloy, J., Leloup, J. C., & Goldbeter, A. (2003). Stochastic models for circadian rhythms: effect of molecular noise on periodic and chaotic behaviour. C. R. Biol., 326 (2), 189–203. Hastings, M. H. (2000). Circadian clockworks: Two loops are better than one. Nat. Rev. Neurosci, 1, 143–146. Khalil, H. K. (1992). Nonlinear Systems. New York, N.Y.: Macmillan. Kholodenko, B. N., Demin, O. V., & Westerhoff, H. V. (1997). Control analysis of periodic phenomena in biological systems. J. Phys. Chem. B., 101 (11), 2070–2081.
D. E. Zak et al. / Computers and Chemical Engineering 29 (2005) 663–673 Larter, R. (1983). Sensitivity analysis of autonomous oscillators: Separation of secular terms and determination of structural stability. J. Phys. Chem., 87 (16), 3114–3121. Larter, R., Rabitz, H., & Kramer, M. (1984). Sensitivity analysis of limit cycles with application to the Brusselator. J. Chem. Phys., 80 (9), 4120– 4128. Leif, F., & Jorgensen, S. B. (2001). Estimation of kinetic parameters in a structured yeast model using regularisation. J. Biotechnol., 88 (3), 223– 237. Leloup, J. C., & Goldbeter, A. (1999). Chaos and birhythmicity in a model for circadian oscillations of the PER and TIM proteins in Drosophila. J. Theor. Biol., 198 (3), 445–459. Leloup, J. C., & Goldbeter, A. (2003). Toward a detailed computational model for the mammalian circadian clock. Proc. Natl. Acad. Sci. USA, 100 (12), 7051–7056. Levsky, J. M., Shenoy, S. M., Pezo, R. C., & Singer, R. H. (2002). Single-cell gene expression profiling. Science, 297 (5582), 836–840. Maly, T., & Petzold, L. R. (1996). Numerical methods and software for sensitivity analysis of differential–algebraic systems. Appl. Numer. Math., 20 (1/2), 57–79. Morohashi, M., Winn, A. E., Borisuk, M. T., Bolouri, H., Doyle, J., & Kitano, H. (2002). Robustness as a measure of plausibility in models of biochemical networks. J. Theor. Biol., 216 (1), 19–30. Novak, B., Pataki, Z., Ciliberto, A., & Tyson, J. J. (2001). Mathematical model of the cell division cycle of fission yeast. Chaos, 11 (1), 277–286. Rabitz, H., & Edelson, D. (1985). Numerical techniques for modelling and analysis of oscillating chemical reactions. In R. J. Field & M. Burger (Eds.), Oscillations and traveling waves in chemical systems (pp.193– 222). John Wiley and Sons. Rao, C. V., Wolf, D. M., & Arkin, A. P. (2002). Control, exploitation and tolerance of intracellular noise. Nature, 420 (6912), 231–237.
673
Savageau, M. A. (1971). Parameter sensitivity as a criterion for evaluating and comparing the performance of biochemical systems. Nature, 229 (5286), 542–544. Schlosser, P. M. (1994). Experimental design for parameter estimation through sensitivity analysis. J. Toxicol. Environ. Health, 43 (4), 495– 530. Shampine, L. F., & Reichelt, M. W. (1997). The MATLAB ODE suite. SIAM J. Sci. Comput., 18 (1), 1–22. Skogestad, S., & Postlethwaite, I. (2000). Multivariable Feedback Control. Chichester, UK: John Wiley and Sons. Stelling, J., & Gilles, E. D. (2001). Robustness vs. identifiability of regulatory modules?. Proc. 2nd Intl. Conf. Systems Biology 181–190. Sveiczer, A., Tyson, J. J., & Novak, B. (2001). A stochastic, molecular model of the fission yeast cell cycle: role of the nucleocytoplasmic ratio in cycle time regulation. Biophys. Chem., 92 (1/2), 1–15. Tomovi´c, R., & Vukobratovi´c, M. (1972). General Sensitivity Theory. New York, N.Y.: American Elsevier. Tyson, J. J., Csikasz-Nagy, A., & Novak, B. (2002). The dynamics of cell cycle regulation. Bioessays, 24 (12), 1095–1109. Vilar, J. M., Kueh, H. Y., Barkai, N., & Leibler, S. (2002). Mechanisms of noise–resistance in genetic oscillators. Proc. Natl. Acad. Sci. USA, 99 (30), 5988–5992. Wang, S. H., & Lee, T. F. (1996). Interpretation of non–minimum phase system with extension to time–varying case. Electron. Lett., 32 (3), 270– 272. Wang, S. H., & Zachery, R. (1996). Singular value decomposition of system input–output matrix and its symmetry property. Comput. Electron. Eng., 22 (3), 231–234. Zak, D. E., Doyle, F. J., Vlachos, D. G., & Schwaber, J. S. (2001). Stochastic kinetic analysis of transcriptional feedback models for circadian rhythms. Proc. 40th IEEE Conf. Decision & Control 849–854.