A method for determining the robustness of biomolecular oscillator models
 Reza Ghaemi^{1},
 Jing Sun^{2},
 Pablo A Iglesias^{3} and
 Domitilla Del Vecchio^{1}Email author
DOI: 10.1186/17520509395
© Ghaemi et al; licensee BioMed Central Ltd. 2009
Received: 13 April 2009
Accepted: 21 September 2009
Published: 21 September 2009
Abstract
Background
Quantifying the robustness of biochemical models is important both for determining the validity of a natural system model and for designing reliable and robust synthetic biochemical networks. Several tools have been proposed in the literature. Unfortunately, multiparameter robustness analysis suffers from computational limitations.
Results
A novel method for quantifying the robustness of oscillatory behavior to parameter perturbations is presented in this paper. This method relies on the combination of Hopf bifurcation and RouthHurwitz stability criterion, which is widely applied in control system design. The proposed method is employed to calculate the robustness of two oscillating biochemical network models previously analyzed in the literature. The robustness bounds here obtained are tighter than what was previously obtained in the literature for both models.
Conclusion
The method here proposed for quantifying the robustness of biochemical oscillator models is computationally less demanding than similar multiparamter variation techniques available in the literature. It also provides tighter bounds on two models previously analyzed in the literature.
Background
It is well acknowledged that characteristics of biochemical systems, such as oscillatory behavior, are preserved under significantly different environmental conditions, i.e., the systems are robust. Some biological systems have been experimentally proven to be robust [1–3]. Since biological systems are robust, the mathematical models developed to represent their characteristics must also reflect this property. In addition, since it is not possible to precisely determine the parameters in the modeling process, biochemical system models must maintain basic properties, such as oscillatory behavior, in the presence of parameter perturbations [4]. Therefore, robustness is considered as an important measure of the validity of mathematical models. Determining the robustness of the oscillatory behavior of a biomolecular model has also relevance in design problems. In fact, by determining the region in parameter space where oscillations persist, one can provide guidelines (see [5, 6] for example) for designing the components of synthetic oscillators such as those of [7, 8].
Parametric robustness of biochemical oscillator networks has been the subject of several studies [4, 9]. In these studies, the maximum allowable parameter deviation from nominal values under which the system oscillations persist is considered as a metric for measuring model robustness. Alternative definitions of parametric robustness are proposed in which the sensitivity of the system equilibrium to parameter variation is considered as a measure of robustness [10]. This measure, however, is only applicable to nonoscillatory systems. Bifurcation analysis is employed to study the sensitivity of oscillatory behavior to variations of a single parameter [11–14]. Systems whose robustness to oneparameteratatime variation is established may be sensitive to simultaneous variation of parameters. Unfortunately, analysis of the effect of multiparameter variation on the oscillatory behavior of a system is more computationally demanding compared to oneparameteratatime variation. In particular, systematic variation of multiple parameters suffers from exponential increase in the number of combinations of parameters to be considered. Therefore, multiparametric sensitivity has often been addressed via computer simulations based on Monte Carlo methods [15]. Since this method relies on random variation of all parameters, the resulting robustness evaluation is inconclusive.
Structured Singular Value (SSV) analysis (μanalysis), a tool developed in the field of robust control, has been employed to provide information on the robustness of systems in the presence of multiple and simultaneous parameter variations from nominal values [11]. In this analysis, a parameter called μ, is calculated whose inverse determines the maximum allowed parameter variation beyond which the system is destabilized. For oscillatory biochemical networks, destabilization means that the system ceases to oscillate. The advantage of μanalysis over Monte Carlo methods is that, due to its deterministic nature, it can compute the extent of parameter uncertainty for which the model is guaranteed to produce the desired behavior. However, for systems with many parameters, μanalysis is not computationally feasible. Hence, one can only rely on computing upper and lower bounds for μ, which determine the maximum stabilizing parameter variation and the minimum destabilizing parameter variation, respectively. To determine whether a system is robust, the lower bound of μ must be calculated. However, available algorithms for computing this lower bound suffer from the curse of dimensionality, i.e., computational time grows exponentially with the number of parameters [9]. Another drawback of μanalysis is that it relies on the linearization of the system about the nominal oscillatory trajectory. As a consequence, a specific combination of parameter variations may destabilize the linearized system while still leading to sustained oscillations in the nonlinear system. This is highly likely in biochemical networks as stable periodic solutions often arise due to nonlinear dynamics.
In addition to μanalysis, analysis of robustness of systems to multiparameter variations has been approached by searching the worst case combination of parameter variations that suppresses oscillations. In [9], the integral of the square of the derivative of one of the states, considered as a measure of the occurrence of oscillations, is minimized with respect to parameter variations employing Hybrid Genetic Algorithms (HGA) [16]. This optimization is performed to determine the region in the parameter space where oscillations persist. Since this method is based on exhaustive simulation, it cannot be applied for systems with large numbers of states and parameters due to the associated prohibitive computational cost. In this paper, we introduce a novel robustness analysis method for oscillatory behavior based on the combination of Hopf bifurcation [17] and RouthHurwitz stability criterion [18]. Combining these two techniques, we compute a scalar parameter ℛ (encompassing all system parameters), which is solely responsible of Hopf bifurcation. As a consequence, we study the persistence of the periodic orbit as this single parameter ℛ is varied. This dramatically reduces the complexity of the problem while retaining the desirable features of multiparameter variation. To translate the maximal variation of ℛ that preserves oscillatory behavior to the maximal variation from a nominal parameter in the original parameter space, we solve an optimization problem. The so obtained maximum parameter variation determines the largest box in parameter space about the nominal parameter values in which the model displays sustained oscillations. Under the assumption that the terms of order higher than three in the Taylor expansion of the system on the center manifold are negligible in the found box, this box provides a tight estimate of the robustness of the system.
We illustrate the application of our approach to models of two molecular networks. The first model describes the molecular network underling adenosine 3',5'cyclic monophosphate (cAMP) oscillations observed in populations of Dictyostelium cells, proposed by Laub and Loomis in [19]. The second model describes the metabolism of an activated neutrophil granulocyte [20].
We next describe the details of the method proposed in this paper.
Methods
in which x(t) ∈ ℝ^{ n }is a vector whose elements are concentrations of chemical species and K ∈ ℝ^{ p }is a vector of parameters. Let K_{0} ∈ ℝ^{ p }be a nominal parameter vector for which system (1) displays a stable periodic orbit. We seek to determine the maximal allowed variation of the parameters K from the nominal value K_{0} before the stable periodic solution disappears. We tackle this problem by analytically computing the set of all values of K at which system (1) admits a Hopf bifurcation. We then compute the largest box in parameter space about K_{0} that does not include any of the values of K at which the system admits a Hopf bifurcation. Under suitable assumptions, this box is the largest box about K_{0} in which sustained oscillations are preserved. This is explained in more detail in the following section.
Hopf bifurcation analysis
Note that even if Hopf bifurcation relies on the properties of the linearization matrix A(K), it allows to infer the existence of a periodic orbit for the original nonlinear system [17].
 (i)
ℛ(K) = 0 if and only if A(K) has two imaginary eigenvalues while all the others having negative real parts;
 (ii)
ℛ (K) > 0 if and only if A(K) has two eigenvalues with positive real parts while all the others having negative real parts;
 (iii)
ℛ(K) < 0 if and only if all eigenvalues of A(K) have negative real parts.
For all K ∈ (K_{0}), ℛ(K) > 0;
There exists on the boundary of B_{δ*}(K_{0}) such that ℛ( ) = 0;
There exists a path Γ ending at in D\B_{δ*}(K_{0}) such that for all points on the path ℛ(K) < 0.
Let K = (k_{1},..., k_{ p }) and assume there exists an open box B_{δ*}(K_{0}) = {K ∈ ℝ^{ p } k_{ i } k_{0, i} <δ* for all i} contained in D such that
According to these assumptions, there is a path Γ in parameter space, crossing through , along which ℛ(K) is first negative, i.e., all eigenvalues of A(K) have negative real parts. Then ℛ(K) = 0 at K = , i.e., A( ) has two imaginary eigenvalues with all the others having negative real parts. Finally, ℛ(K) > 0 in B_{δ*}(K_{0}), i.e., A(K) has two eigenvalues with positive real parts with all the others having negative real parts. Therefore, a Hopf bifurcation occurs at , corresponding to ℛ( ) = 0. A small amplitude periodic orbit appears for small ℛ(K) > 0 if the Hopf bifurcation is supercritical. Hence, if ℛ(K) is sufficiently small for all K ∈ B_{δ*}(K_{0}) the system admits a stable periodic orbit in B_{δ*}(K_{0}) (having ℛ(K) small enough in the box B_{δ*}(K_{0}) implies that the Taylor expansion of the system on the center manifold passing through and x_{ e }( ) is such that the terms of order higher than three are negligible in B_{δ*}(K_{0})). According to Hopf bifurcation theorem, at K = , the periodic orbit disappears and therefore, we take d* as a robustness measure.
The assumption that ℛ(K) is sufficiently small in the open box δ* (K_{0}) may not hold for all biochemical systems. One can further validate that ℛ(K) is sufficiently small in the open box B_{δ*}(K_{0}) by simulating the system for ℛ(K) between zero and its maximum in B_{δ*}(K_{0}) and by verifying that oscillations persist. Since simulation is performed as one scalar parameter ℛ(K) varies (as opposed to varying multiple parameters at once), it is computationally feasible. It is worth noting that there may be variations in how the higher order terms in the Taylor expansion grow when ℛ(K) grows depending on the path adopted to vary ℛ(K) inside the box B_{δ*}(K_{0}).
Determining the Rfunction ℛ(K)
Tabular method for a system with characteristic polynominal C (s, k)
s ^{ n }  a _{0}  a _{2}  a _{4}  ⋯ 
s ^{n1}  a _{1}  a _{3}  a _{5}  ⋯ 
s ^{n2} 


 ⋯ 
s ^{n3} 


 ⋯ 
⋮  ⋮  ⋮  ⋮  ⋱ 
s ^{3} 

 0  ⋯ 
s ^{2} 
 0  ⋯  
s ^{1} 
 0  0  ⋯ 
s ^{0} 
 0  0  ⋯ 
as the Rfunction.
Algorithm for calculating the maximum allowed uncertainty δ*
In light of the earlier section, the resulting algorithm for calculating the maximum allowed uncertainty δ* about K_{0} for which the stable periodic orbit is preserved is summarized as follows.
Algorithm 1
Step 1. Calculate the maximum value of δ*, such that for all K ∈ B_{δ*}(K_{0}), we have that ℛ(K) > 0.
Step 2. Identify on the boundary of B_{δ*}(K_{0}) such that ℛ( ) = 0.
Step 3. Verify that a_{0}(K), a_{1}(K), (K, ⋯, (K), and (K), computed in equation (4), are positive for all K ∈ B_{δ* + ϵ}(K_{0}) for some ϵ > 0.
Step 4. Verify that in an open neighborhood N( ) ⊂ B_{δ* + ϵ}(K_{0}) about , we have that ℛ(K) ≤ 0 implies (K) is positive. Set D := N( ) ∩ B_{δ*}(K_{0}).
Step 5. Verify that the value of ℛ(K) for K on the path Γ = {(1 α)K_{0} + α  α > 1} ∩ B_{δ* + ϵ}(K_{0}) is negative.
Remark 1. If a_{0}(K), a_{1}(K), (K, ⋯, (K), and (K) are all positive in D, Algorithm 1 can be employed with ℛ(K) :=  (K).
This optimization problem provides the value of δ* of Step 1 and the value of of Step 2. Similarly, Steps 3 and 4 can be implemented by determining the minimum of the designated functions on the specified areas in the parameter space. Step 5 can be implemented through a scalar optimization problem, which can be solved via the bisection method [23].
Remark 2. The robustness of the periodic orbit of a system can be evaluated by determining the deviation δ* from the nominal parameter values at which the function ℛ changes sign. This deviation can be numerically computed by employing standard optimization techniques such as SQP and HGA. As a consequence, this technique is computationally lighter than multiparameter robustness analysis based on randomsearch methods, in which the system is simulated at each point in parameter space which is determined by the random search algorithm.
Remark 3. If some key parameters are a priori identified to have more significant influence on the system dynamics than others, the proposed algorithm can be employed to analyze the effect of variation of significant parameters while nonsignificant ones are kept constant at their nominal values. A method for determining the robustness interval specifically for individual parameters is considered in [9].
Results
In this section, we illustrate the detailed application of Algorithm 1 to two well known models of oscillatory biochemical networks and compare the results with robustness measures previously obtained in the literature.
Application to the Laub and Loomis model
Nominal values for each parameter
Parameter  Units  Nominal Value 

k _{0,1}  min ^{1}  2.0 
k _{0,2}  Mol ^{1} min ^{1}  0.9 
k _{0,3}  min ^{1}  2.5 
k _{0,4}  min ^{1}  1.5 
k _{0,5}  min ^{1}  0.6 
k _{0,6}  Mol ^{1} min ^{1}  0.8 
k _{0,7}  Mol ^{1} min ^{1}  1.0 
k _{0,8}  Mol ^{1} min ^{1}  1.3 
k _{0,9}  min ^{1}  0.3 
k _{0,10}  Mol ^{1} min ^{1}  0.8 
k _{0,11}  min ^{1}  0.7 
k _{0,12}  min ^{1}  4.9 
k _{0,13}  min ^{1}  23.0 
k _{0,14}  min ^{1}  4.5 
Equilibria of the system
For the detailed derivation, the reader is referred to Appendix 2.
Determining the function ℛ
Routh Hurwitz table for characteristic polynomial (9)
s^{7}  1  a_{2}  a_{4}  a_{6} 
s ^{6}  a_{1}  a_{3}  a_{5}  a_{7} 
s ^{5} 


 0 
s ^{4} 


 0 
s ^{3} 

 0  0 
s ^{2} 

 0  0 
s ^{1} 
 0  0  0 
s ^{0} 
 0  0  0 
Implementation of Algorithm 1
where Δ = (δ_{1}, ⋯, δ_{14}). This problem is solved numerically employing the Sequential Quadratic Programming (SQP) method [24]. Solving the optimization problem (12) results in δ* = 0.0051 = Δ*∞, where . Although the SQP provides only a local optimal solution, since δ* is a small number, we expect it to be the actual optimal solution. To further verify that the solution is the global optimal, we perform exhaustive search using gridding methods for the points {k_{0, i} δ*k_{0, i}, k_{0, i}, k_{0, i}, + δ*k_{0, i}}. To perform exhaustive search we evaluate ℛ(K) for 3^{14} values of K and determine the minimum value. The resulting minimum value is the same as the one obtained with the SQP.
Elements of the vector
Parameter  Units  Perturbed Value 

 min ^{1}  1.9898 
 Mol ^{1} min ^{1}  0.8954 
 min ^{1}  2.5128 
 min ^{1}  1.5076 
 min ^{1}  0.5969 
 Mol ^{1} min ^{1}  0.8041 
 Mol ^{1} min ^{1}  1.0051 
 Mol ^{1} min ^{1}  1.2934 
 min ^{1}  0.3015 
 Mol ^{1} min ^{1}  0.8041 
 min ^{1}  0.6964 
 min ^{1}  4.9250 
 min ^{1}  22.8827 
 min ^{1}  4.5229 
where Δ = (δ_{1}, ⋯, δ_{14}).
The optimization problem is solved for j = 0, ⋯, 6 using the gradient descent method. We obtain the minimum values 15.38, 76.13, 184.64, 225.73, 128.74, 93.87 for the functions and , respectively. Another approach for calculating the minimum value of the functions is performing exhaustive search in the gridded space and we consider the elements
{k_{0, i} (δ* + ϵ)δk_{0, i}, k_{0, i}, k_{0, i}+ (δ* + ϵ)k_{0, i}i = 1, ⋯, 14}. Therefore, to perform exhaustive search we evaluate each of the functions , j = 0, ⋯, 6 for 3^{14} elements (3 elements for each dimension) and determine the minimum value of each function. The resulting minimum values are the same as those obtained employing the gradient descent method.
Since in the box B_{δ* + ϵ}(K_{0}) all elements of the vector v are positive, the conditions described in Steps 3 and 4 are verified in D = B_{δ* + ϵ}(K_{0}), according to Remark 1.
Step 5. Evaluating ℛ(K) for K on the path {(1  α)K_{0} + α  α > 0} ∩ B_{δ* + ϵ}(K_{0}), we can show that it is negative except for at which it is zero and Step 5 of Algorithm 1 is completed.
Therefore, all the Steps in Algorithm 1 are completed and the system has a stable periodic orbit in the open box B_{δ*}(K_{0}) with δ* = 0.0051.
To further verify that ℛ is sufficiently small in B_{δ*}(K_{0}), we simulate the system for ℛ between zero and its maximum in B_{δ*}(K_{0}) and verify that oscillations persist (see Appendix 3 for the details).
Application to the model of the metabolism of activated neutrophils
Reaction and rate constants of neutrophils
Reaction  Rate expression (R_{ i })  Rate constant 

R_{1}(x)  k_{1}x_{5}x_{1}  k_{1}x_{2}  K_{1} = 5.0 × 10^{7}M^{}1s^{1} 
k_{1} = 58s^{1}  
R_{2}(x)  k _{2} x _{2} x _{8}  k_{2} = 1.0 × 10^{7}M^{}1s^{1} 
R_{3}(x)  k _{3} x _{3} x _{8}  k_{3} = 4.0 × 10^{3}M^{}1s^{1} 
R_{4}(x)  k _{4} x _{1} x _{6}  k_{4} = 2.0 × 10^{7}M^{}1s^{1} 
R_{5}(x) 
 k_{5} = 1.0 × 10^{7}M^{}1s^{1} 
R_{6}(x)  k _{6} x _{4} x _{6}  k_{6} = 1.0 × 10^{5}M^{}1s^{1} 
R_{7}(x)  k _{7} x _{10} x _{14}  k_{7} = 1M^{}1s^{1} 
R_{8}(x)  k _{8} x _{11} x _{14}  k_{8} = 5.0 × 10^{7}M^{}1s^{1} 
R_{9}(x) 
 k_{9} = 5.0 × 10^{8}M^{}1s^{1} 
R_{10}(x)  k _{10} x _{16} x _{10}  k_{10} = 1.0 × 10^{7}M^{}1s^{1} 
R_{11}(x) 
 k_{11} = 6.0 × 10^{7}M^{}1s^{1} 
R_{12}(x)  k _{12}  k_{12} = 3.0 × 10^{5}M^{}1s^{1} 
R_{13}(x)  k_{13} k_{13x 14}  k_{13} = 1.25 × 10^{5}M^{}1s^{1} 
k_{13} = 4.5 × 10^{2}M^{}1s^{1}  
R_{14}(x)  k_{14}(x_{7}x_{14})  k_{14} = 30s^{1} 
R_{15}(x)  k_{15}(x_{5}x_{12})  k_{15} = 30s^{1} 
R_{16}(x)  k_{16}(x_{8}x_{15})  k_{16} = 10s^{1} 
R_{17}(x)  k_{17}(x_{9}x_{16})  k_{17} = 10s^{1} 
R_{18}(x)  k_{18}(x_{6}x_{13})  k_{18} = 0.01s^{1} 
R_{19}(x) 
 V = 288 × 10^{6}M^{}1s^{1} 
L = 550  
k_{0} = 1.5 × 10^{6}M  
k_{ N }= 60 × 10^{6}M 
Equilibria of the system
The system model does not belong to the category of Ssystems. Therefore, we use Newton's method to calculate the equilibrium, x_{ e }(K), of system (15) by solving f(x_{ e }, K) = 0 for x_{ e }. We first simulate the system with nominal parameters K_{0} and choose a point on the corresponding periodic orbit to initialize the Newton iterations to calculate x_{ e }(K_{0}). Then for any parameter vector K ≠ K_{0}, the Newton's method is initialized by x_{ e }(K_{0}) and the iteration is performed to achieve x_{ e }(K).
Determining the function ℛ
Routh Hurwitz table
s ^{14}  1  a _{2}  a _{4}  ⋯  a _{14} 
s ^{13}  a _{1}  a _{3}  a _{5}  ⋯  0 
s ^{12} 


 ⋯  
⋮  ⋮  ⋮  ⋮  ⋱  ⋮ 
s ^{3} 

 0  ⋯  
s ^{2} 

 0  ⋯  
s ^{1} 
 0  0  ⋯  0 
s ^{0} 
 0  0  ⋯  0 
Implementation of Algorithm 1
We now apply Algorithm 1 to determine the region in which the stable periodic orbit persists.
This problem is equivalent to solving the following optimization problem: δ subject to: δ_{ i } ≤ δ, i = 1, ⋯, 26 and ℛ(k_{0,1} + δ_{1}k_{0,1}, ⋯, k_{0,26} + δ_{26}k_{0,26}) ≤ 0. Employing the SQP solver, we achieve δ* ≃ .17. However, employing first Adaptive Search Algorithm [27] and then SQP to the result of the Adaptive Search Algorithm (see Appendix 4 for the details), we achieve δ* = .0591. HGA and then SQP are also employed to calculate δ* which leads to the same results as the algorithm described in Appendix 4.
Step 2. The solution to the above optimization problem provides the optimal values for δ_{ i }, i = 1, ⋯, 14, i.e., that make the constraint ℛ(K) ≥ 0 active. Therefore, at the point defined as in which [δ_{1}, ⋯, δ_{26}] = 0.0591 [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0], which is on the boundary of the box B_{δ*}(K_{0}), ℛ vanishes, i.e., ℛ( ) = 0. It can be shown that at the point , with ϵ being sufficiently small, the linearized system about the equilibrium is stable and the system ceases to oscillate.
Step 3. Let us consider the box B_{δ* + ϵ}(K_{0}) for ϵ = 0.0001. Then, in B_{δ* + ϵ}(K_{0}) the functions and are positive if their minimum in the box is positive. We determine the minimum of each of the functions in the box B_{δ* + ϵ}(K_{0}) employing a HGA and SQP. The minimum value of the functions and are found to be 172.91, 8214.8, 173070, 1.8478e + 006, 1.0279e + 007, 2.7101e + 007, 3.6251e + 007, 2.6068e + 007, 1.0073e + 007, 1.7745e + 006, 1.3174e + 004, and 5.1522e  004, respectively.
Step 4. By defining a small neighborhood about , that is, B_{.001}( ), we obtain that (K) is positive in B_{.001}( ) and therefore Step 4 of Algorithm 1 is completed with N( ) := B_{.001}( ) and D := N( ) ∪ B_{δ*}(K_{0}).
Step 5. We evaluate ℛ(K) for K on the path Γ = {(1  αK_{0}) + α  α > 0} ∩ B_{δ* + ϵ}(K_{0}) and confirm that it is negative except on where ℛ is zero. Therefore, Step 5 of Algorithm 1 is completed and the largest box about K_{0} in which the system has a stable periodic orbit is B_{δ*}(K_{0}) with δ* = 0.0591, provided that ℛ is sufficiently small on B_{δ*}(K_{0}).
To verify that ℛ is sufficiently small in the box B_{δ*}(K_{0}), we simulate the system for ℛ ranging from zero to the maximum attained in B_{δ*}(K_{0}) and verify that oscillations persist (see Appendix 5 for the details).
Discussion
The method introduced in this paper relies on the computation of ℛ, that is, the scalar function of system parameters whose sign determines the existence of the stable periodic orbit. The robustness of the periodic orbit of a system can be evaluated by determining the deviation δ* from the nominal parameter values at which the function ℛ changes sign. This deviation can be numerically computed by employing standard optimization techniques such as SQP and HGA. As a consequence, this technique is computationally lighter than multiparameter robustness analysis based on randomsearch methods, in which the system is simulated at each point in parameter space. Nevertheless, since ℛ encompasses all system parameters, this method retains the desirable features of multiparameter robustness analysis. Moreover, the proposed method provides more accurate robustness measures when compared to methods based on the linear approximation of the system about the nominal periodic orbit, as it is shown in the first application example. The proposed method relies on Hopf bifurcation theorem and on the assumption that ℛ(K) is sufficiently small for K inside the box B_{δ*}(K_{0}). This guarantees that the terms of order higher than three in the Taylor expansion of the system on the center manifold about and x_{ e }( ) are negligible in the box B_{δ*}(K_{0}). To provide evidence that this assumption is satisfied, we vary ℛ from its minimum to its maximum in the box B_{δ*}(K_{0}) and verify via simulation that the periodic orbit persists. This simulation step does not present computational limitations. It is in fact performed as one scalar parameter (ℛ) is varied as opposed to varying multiple parameters at once as in multiparameter robustness analysis. The underlying assumption for the approach presented in this paper to provide a tight bound on the robustness of the system is that the nominal periodic orbit of the system originates from a Hopf bifurcation. If this were not the case, the provided bound would not necessarily be meaningful as other types of bifurcations may be responsible of the birth and death of the periodic orbit. Therefore, the approach of this paper is generally applicable and restricted to those natural oscillatory systems exhibiting Hopf bifurcation.
For the Laub and Loomis model, previous work employing HGA, in which the system is simulated at each point in parameter space, the robustness of the system was determined as 0.6% [9]. Employing the method proposed in this paper, the robustness of the system has been determined as δ* = 0.51%. This bound is tight, as the system ceases to oscillate at a combination of parameters that is away from the nominal value K_{0} only slightly more than δ*. This bound is therefore tighter than the one found employing μ analysis or global/hybrid optimization methods.
For the model of the oscillatory metabolism of activated neutrophils, previous work only performed oneparameteratatime variation [12]. According to this single parameter variation analysis, the minimum deviation from nominal values of parameters which causes the periodic orbit to disappear is 16.67%. By contrast, employing the method proposed in this paper, the robustness of the system has been quantified as δ* = 5.91%. This bound is tight as a combination of parameters that is away from the nominal point K_{0} slightly more than δ* has been determined at which the system ceases to oscillate. This result shows that the oscillatory behavior of this model is not as robust with respect to parameter variation as it was perceived.
Conclusion
The robustness analysis of biomolecular systems is an important problem in systems and synthetic biology. Previously, the robustness of a system with respect to parameter variations has been investigated by employing μ analysis on the linearized system about nominal periodic orbit or by applying HGA, in which the system is simulated for each combination of parameter values. In this paper, a method based on the combined application of Hopf bifurcation and RouthHurwitz stability criterion is introduced. We computed a scalar function of all system parameters whose sign determines the existence of a stable periodic orbit. This method is applied to two biomolecular systems: the Laub and Loomis model and the model of the oscillatory metabolism of activated neutrophils. The maximum allowed parameter variation with respect to nominal values under which the system preserves oscillatory behavior is calculated. For the Laub and Loomis model, the computed maximum allowed variation is tighter than what was obtained with previous multiparametric analysis methods. For the model of activated neutrophils, only single parameter variations were considered in the literature to evaluate parametric robustness. Employing the method proposed in this paper, we evaluate the robustness of the system to be about one third of the one estimated in the literature employing single parameter variations.
Appendix
Appendix 1
The following result gives conditions for which T(K) can be taken as an Rfunction.
Claim 1. Assume that for all K ∈ D,
(a) a_{0}(K), a_{1}(K), (K), ⋯, (K) and (K) are positive;
(b) K ∈ D ∩ {K  T(K) ≥ 0} implies that (K) > 0.
Then ℛ(K) = T(K) is an Rfunction.
Proof. If T(K) = (K) (K) < 0 for K ∈ D, then when K ∈ B_{δ*}(K_{0}) either (K) is positive and (K) is negative or viceversa. Therefore, given (a), there are always two changes of sign in the vector v(K) in (6) for all K ∈ D. Hence, according to RouthHurwitz criterion, there are two eigenvalues with positive real part while all the others have negative real part in D. Therefore, ℛ(K) = T(K) satisfies property (ii). According to condition (b), T(K) = 0 implies (K) > 0 and (K) = 0. Therefore, from the RouthHurwitz criterion, there are two imaginary eigenvalues and all the others have negative real part. As a consequence, ℛ(K) = T(K) satisfies property (i).
According to (b), if T(K) > 0 when K ∈ D, then (K) > 0. This, along with condition (a), implies by Routh Hurwitz criterion that all eigenvalues of A(K) have negative real parts. Consequently, ℛ(K) = T(K) satisfies property (iii).
Appendix 2
Therefore, the system has two equilibrium points. One is determined by equation (25), which we do not consider for Hopf bifurcation, because negative values for the states x_{1}, ⋯, x_{7} would result from oscillations and protein concentration is always positive. The other nonzero equilibrium, employed next, is determined by equations (24262930313233).
Appendix 3
In order to perform this onedimensional parameter study, we define a path in parameter space along which ℛ changes from zero to its maximum value in the box B_{δ*}(K_{0}) monotonically. In the following, we show that ℛ is monotonically increasing with respect to parameter k_{1} in B_{δ*}(K_{0}). Therefore, the path is created by varying k_{1} from zero at to the point at which ℛ reaches its maximum value, while keeping all other parameters constant to the values taken at .
in which Δ = [δ_{1}, ⋯, δ_{14}]. Employing the gradient descent method, we obtain the maximum value of H(K) over the box B_{δ*}(K_{0}) as 24.9152. Therefore, ℛ(K) is positive over the box B_{δ*}(K_{0}), which implies that ℛ is monotonically increasing with respect to k_{1}.
where Δ = [δ_{1}, ⋯, δ_{14}]. This problem is solved again through the gradient decent method. Denoting the solution of (36), we have that = 4.9828. By varying parameter k_{1}, we thus simulate the system for ℛ ∈ [0, ]. Simulation results show that the system has a periodic orbit for 0 < ℛ≤ .
Appendix 4
Consider the optimization problem δ subject to: δ_{ i } ≤ δ, i = 1, ⋯, 26 and ℛ(k_{0,1} + δ_{1}k_{0,1}, ⋯, k_{0,26} + δ_{26}k_{0,26}) ≤ 0. We adopt the following solution:
1 Set δ^{ l }= 0.17, l = 1, K_{ init }= K_{0}.
We solve this optimization problem using the direct search method of the Genetic Algorithm toolbox of Matlab, with mesh adaptive search algorithm [27], using the initial guess K_{ init }. The solution is taken as an initial guess for SQP and the solution provided by SQP is considered as the solution of the above optimization problem. If at the optimal point , ℛ(K^{ l }) is nonnegative, terminate.
3 Determine 0 <α ≤ 1 such that ℛ(αK_{0} + (1  α)K^{ l }) = 0 (using methods such as bisection [23]).
4 Set K_{ init }= αK_{0} + (1 α)K^{ l }, l = l + 1 and .
5 Go to 2.
Appendix 5
In order to simulate the system over different values of ℛ, we vary the value of α from 0 to 1 at 855 points such that ℛvaries of at most 20 at each step. The simulation results show that the system has a periodic orbit over the interval 0 < ℛ < 18665.
Declarations
Acknowledgements
This work was in part supported by The Rackham School at University of Michigan.
Authors’ Affiliations
References
 Barkai N, Leibler S: "Robustness in simple biochemical networks". Nature. 1997, 387: 913917. 10.1038/43199View ArticlePubMedGoogle Scholar
 Alon U, Surette M, Barkai N, Leibler S: "Robustness in bacterial chemotaxis". Nature. 1999, 397: 168171. 10.1038/16483View ArticlePubMedGoogle Scholar
 von Dassow G, Meir E, Munro EM, Odell GM: "The segment polarity network is a robust development module". Nature. 2000, 406: 188192. 10.1038/35018085View ArticlePubMedGoogle Scholar
 Savageau MA: " Parameter sensitivity as a criterion for evaluating and comparing the performance of biochemical systems". Nature. 1971, 229: 542544. 10.1038/229542a0View ArticlePubMedGoogle Scholar
 Del Vecchio D: "Design and Analysis of an ActivatorRepressor Clock in E. coli". Proceedings of American Control Conference, New York. 2007Google Scholar
 El Samad H, Del Vecchio D, Khammash M: "Repressilators and Promotilators: Loop Dynamics in Gene Regulatory Networks". Proceedings of American Control Conference, Santa Barbara. 2005Google Scholar
 Elowitz MB, Leibler S: "A Synthetic Oscillatory Network of Transcriptional Regulators". Nature. 2000, 403: 335338. 10.1038/35002125View ArticlePubMedGoogle Scholar
 Atkinson M, Savageau M, Myers J, Ninfa A: "Development of Genetic Circuitry Exhibiting Toggle Switch or Oscillatory Behavior in Escherichia coli". Cell. 113 (5): 597607.
 Kim J, Bates DG, Postlethwaite I, Ma L, Iglesias PA: "Robustness analysis of biochemical network models". IEE Proceedings Systems Biology. 2006, 153 (3): 96104. 10.1049/ipsyb:20050024View ArticlePubMedGoogle Scholar
 Chen BS, Wang YC, Wu WS, Li WH: "A new measure of the robustness of biochemical networks". Bioinformatics. 2005, 21 (11): 26982705. 10.1093/bioinformatics/bti348View ArticlePubMedGoogle Scholar
 Ma L, Iglesias PA: "Quantifying robustnes of biochemical network models". BMC Bioinformatics. 2002, 3: 3850. 10.1186/14712105338PubMed CentralView ArticlePubMedGoogle Scholar
 Jacobson EW, Cedersund G: "Structural robustness of biochemical network modelswith application to the oscillatory metabolism of activated neutrophils". IET Systems Biology. 2008, 2 (1): 3947. 10.1049/ietsyb:20070008View ArticleGoogle Scholar
 Wolf J, BeckerWeiman S, Heinrich R: "Analysing the robustness of cellular rhythms". IEE Proceedings Systems Biology. 2005, 2 (1): 3541. 10.1049/sb:20045035View ArticlePubMedGoogle Scholar
 Stelling J, Gilles ED, Doyle FJ: "Robustness properties of circadian clock architectures". Proc Natl Acad Sci USA. 3215, 101 (36): 132101. 10.1073/pnas.0401463101.View ArticleGoogle Scholar
 Liu JS: "Monte Carlo strategies in scientific computing". 2001, New York: SpringerGoogle Scholar
 Lobo FG, Goldberg DL: "Decision making in a hybrid genetic algorithm". 1996, Technical Report No. 96009, Illinois Genetic Algorithms LaboratoryGoogle Scholar
 Wiggins S: "Introduction to Applied Nonlinear Dynamical Systems and Chaos". 2003, New York: SpringerGoogle Scholar
 Hurwitz A: "On the conditions under which an equation has only roots with negative real parts". Selected Papers on Mathematical Trends in Control Theory, New York, Dover. 1964Google Scholar
 Laub MT, Loomis WF: "A molecular network that produces spantaneous oscillations in excitable cells of Dictyostelium". Mol Biol Cell. 1998, 9: 35213532.PubMed CentralView ArticlePubMedGoogle Scholar
 Olsen LF, Kummer U, Kindzelskii AL, Petty HR: "A model of oscillatory metabolism of activated neutrophils". Biophysical Journal. 2003, 84 (4): 6981. 10.1016/S00063495(03)748334PubMed CentralView ArticlePubMedGoogle Scholar
 Voit EO: "Computational Analysis of Biochemical System". 2000, Cambridge University Press, Cambridge, UKGoogle Scholar
 Kelley CT: "Solving Nonlinear Equations with Newton's Method". Fundamentals of Algorithms, SIAM, Philadelphia. 2003Google Scholar
 Burden RL, Faires JD: "Numerical Analysis". 2000, Thomson Brooks/cole; Florence, KYGoogle Scholar
 Flecher R: "Practical methods of Optimization". 1981, 2: Constrained Optimization, Wiley, Chichester, EnglandGoogle Scholar
 Petty RH: "Neutrophil oscillations: temporal and spatiotemporal aspects of cell behaviour". Immunol Res. 2000, 23: 8594. 10.1385/IR:23:1:85.View ArticleGoogle Scholar
 Amit A, Kindzelskii AL, Zanont J, Jarvis JN, Petty RH: "Complement deposition of immune complexes reduces the frequencies of metabolic, proteolytic and superoxide oscillations of migrating neutrophils". Cell Immunol. 1999, 194 (1): 4753. 10.1006/cimm.1999.1481View ArticlePubMedGoogle Scholar
 Michael LR, Torczon V: "Pattern Search Methods for Linearly Constrained Minimization". SIAM Journal on Optimization. 2000, 10 (3): 916941.Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.