 Research article
 Open Access
 Published:
A method for determining the robustness of biomolecular oscillator models
BMC Systems Biology volume 3, Article number: 95 (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 this section, we propose a robustness analysis technique based on the combination of Hopf bifurcation [17] and RouthHurwitz stability criterion [18]. Consider a model of a biochemical network given as
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
At a Hopf bifurcation, an equilibrium of a dynamical system loses stability as a pair of complex conjugate eigenvalues of the linearization about the equilibrium cross the imaginary axis of the complex plane. We thus first compute the equilibrium of the system as a function of the parameters, i.e., x_{ e }(K) such that f(x_{ e }(K), K) = 0. In systems with Ssystem representation [10, 21], the equilibrium can be calculated analytically. For general systems, we rely on numerical methods such as Newton's iterations to compute x_{ e }(K) [22]. Once a nominal equilibrium is identified, i.e., x_{ e }(K_{0}), we initialize Newton's method with this nominal equilibrium to calculate the equilibrium when K ≠ K_{0}. Then, to determine conditions under which the system undergoes a Hopf bifurcation, we study the behavior of the eigenvalues of the matrix
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].
Assume there exists an open region D ⊆ ℝ^{p}in parameter space about the nominal parameter vector K_{0} = (k_{0, 1},..., k_{0, p}) and a function ℛ: ℝ^{p}→ ℝ, called Rfunction, with the following properties:

(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
The sets D, B_{δ*}(K_{0}), and the path Γ are illustrated in Figure 1.
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)
To construct the function ℛ, we consider the characteristic polynomial of A(K)
in which the coefficients of the polynomial are functions of the parameter vector K. To evaluate the effect of K on the eigenvalues of A(K), we appeal to the well known RouthHurwitz stability criterion [18]. This criterion can be translated into a tabular method, in which for a system with characteristic polynomial C(s, K), the table has n + 1 rows and the following structure depicted in Table 1 (neglecting the dependence on K). In this table, := a_{2k1}, := a_{2(k1)}, k = 1, 2, ⋯, ,
and
According to the RouthHurwitz criterion, the number of eigenvalues of A(K) with positive real part is determined by the number of sign changes in the vector
As shown in Appendix 1, we can take
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).
Step 5 verifies that ℛ is negative on the path Γ outside the region B_{δ*}(K_{0}). Steps 3 and 4 verify that the conditions under which ℛ(K) = T(K) as given in Appendix 1 are satisfied. In particular, since ℛ(K) ≤ 0 is equivalent to T(K) ≥ 0, Step 4 verifies that condition (b) of Claim 1 in Appendix 1 is satisfied. Step 3 verifies that condition (a) of Claim 1 in Appendix 1 is satisfied. Once the function ℛ(K) has been computed, Step 1 can be implemented by calculating the smallest δ for which ℛ(K) is nonpositive on some point in the closed box B_{ δ }(K_{0}). This can be obtained by solving the following optimization problem
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
First, we consider a model of the molecular network underlying adenosine 3',5'cyclic monophosphate (cAMP) oscillations observed in populations of Dictyostelium cells, proposed by Laub and Loomis in [19]. The model, based on the network depicted in Figure 2, displays spontaneous oscillations in cAMP observed during the early development of Dictyostelium discoideum.
In this model, changes in the enzymatic activities of these proteins are described by the following system of seven nonlinear differential equations
in which the state variable x = [x_{1}, ⋯, x_{7}] represents the concentration of the seven proteins: x_{1} = [ACA], x_{2} = [PKA], x_{3} = [ERK2], x_{4} = [REGA], x_{5} = [Internal cAMP], x_{6} = [External cAMP] and x_{7} = [CAR1]. The fourteen coefficients, k_{ i }, i = 1, ⋯, 14, are system parameters and we denote K = (k_{1}, ⋯, k_{14}). It is shown in [19] that oscillations appear at the nominal parameter values given in Table 2. We define the nominal parameter vector K_{0} := [k_{0, 1}, ⋯, k_{0, 14}]. To calculate the matrix A(K) in equation (2), we compute in the next section the equilibrium of system (7).
Equilibria of the system
Since the system has Ssystem structure [10], i.e., the elements of f(x, K) are the addition of two monomials, the equilibrium can be calculated analytically. The equilibrium we consider is given by , in which,
For the detailed derivation, the reader is referred to Appendix 2.
Determining the function ℛ
In this section, we employ the explicit representation of the equilibria computed in the previous section to linearize the system and compute the function ℛ. The linearization of the system about the nonzero equilibrium renders the linearization matrix
The characteristic polynomial of the matrix A(K) can be written in the following form (neglecting the dependence on K)
in which a_{ i }, i = 0, ⋯, 6, are calculated analytically as functions of k_{1}, ⋯, k_{14}. The corresponding RouthHurwitz table is given in Table 3.
Following Remark 1, we set ℛ(K) = and then show that the values of a_{1}(K), are all positive in the found set D. Since the expression of the function ℛ(K) is very long, its exact formula is omitted here. Instead, to provide a qualitative understanding of the behavior of ℛ as a function of the parameter vector K, Figure 3 shows the 2dimensional zero level set of the function ℛ for different nominal values of K.
Implementation of Algorithm 1
Step1. To implement Step 1 of Algorithm 1, we determine the largest δ such that ℛ(K) is positive in the open box B_{ δ }(K_{0}). The box is defined as follows to reflect the relative variation of parameters instead of absolute change:
We determine the maximum percentage of the parameter uncertainty under which ℛ(K) ≥ 0 via solving the following optimization problem:
in which Δ = (δ_{1}, ⋯, δ_{14}) and k_{0,1}, ⋯, k_{0,14} are the nominal parameters. The above optimization problem determines the smallest δ such that the box B_{ δ }(K_{0}) contains an element at which ℛ is nonpositive. This is equivalent to determining the largest δ such that ℛ is positive inside B_{ δ }(K_{0}). The optimization problem (11) has the following equivalent form:
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.
Step 2. The solution to the optimization problem (12) provides the optimal values of δ_{ i }, i = 1, ⋯, 14, i.e., , which make the constraint ℛ(K) ≤ 0 active. Therefore, at the point defined as
which is on the boundary of the box B_{δ*}(K_{0}), we have that ℛ() = 0 and at least for some 1 ≤ i ≤ 14  = δ*. The numerical values of the vector are provided in Table 4. It can be shown that at , with ϵ being sufficiently small, the linearized system around the equilibrium is stable and the system ceases to oscillate.
Steps 3 and 4. To implement Steps 3 and 4 of Algorithm 1, we consider Remark 1 and evaluate the sign of the elements in the vector for all K in the set D defined as D := B_{δ* + ϵ}(K_{0}), where ϵ = 0.0001. The sign of these parameters are evaluated via solving the following optimization problem for each function in . This leads to the minimum value that these functions take in B_{δ* + ϵ}(K_{0}). Defining := a_{1}, this minimum value is provided by
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
The second model we consider for robustness analysis describes the metabolism of an activated neutrophil granulocyte. Neutrophils constitute the pivotal part of the defence system against invading pathogens or bacteria. Upon bacterial invasion, neutrophils leave bloodstream and migrate actively toward the site of infection where they absorb and kill the bacteria. The necessary antibacterial and digestive materials are produced upon activation of the neutrophil. The activation dramatically increases the production of the reduced form of nicotinamide adenine dinucleotide phosphate (NADPH) via hexose monophosphate shunt and initiates the production of the NADPH oxidase complex that assembles at the phagosomal membrane. Electrons are transferred from cytoplasmic NADPH to oxygen on the phagosomal side of the membrane, generating the socalled reactive oxygen species by creating superoxide as an intermediate step [20, 12]. It is shown that in migrating neutrophils the concentration of NADPH and reactive oxygen species oscillate [25, 26]. In this paper, we consider the model of the system that is presented by Olsen et al. [20]. The model has 16 states and 24 parameters and it can be written as
in which the state variable x = [x_{1}, ⋯, x_{16}] represents the concentration of the 16 species. The factor ρ refers to the fractional volume of the phagosome over the cytosol and is assumed to be 0.1. The functions R_{1},..., R_{19} represent reaction rates and are described in Table 5 along with nominal parameter values.
In model (15), there are only 14 independent states as and . Therefore, we let
Substituting x_{1} and x_{16} in (15) and employing (16), we have (with abuse of notation) the new state vector x = [x_{2}, ⋯, x_{15}] with independent states and new differential equation
in which the parameter vector K ∈ ℝ^{26} includes all 24 parameters in Table 5 in addition to k_{25} and k_{26}, which depend on the initial concentration of the molecular species. Assuming initial concentration of 300 μM for ferric peroxidase (x_{1}) and melatonin (x_{16}) and zero for the rest, we have a nominal value of 300 μ M for k_{25} and k_{26} as in [20]. Moreover, the nominal parameter vector is defined as
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 ℛ
The linearization of the system about the nonzero equilibrium renders the matrix A ∈ ℝ^{14 × 14}, which is calculated numerically as a function of the equilibrium x_{ e }and parameter vector K. The 14th order characteristic polynomial of the matrix A is given (omitting the dependence on K) by C(s) = s^{14} + a_{1}s^{13} + ⋯ + a_{13}s + a_{14}. Table 6 shows the RouthHurwitz table in which the parameters are functions of x_{ e }and K. Using Newton's method, x_{ e }is calculated as a function of K and therefore the vector is also calculated as a function of K. According to equations (5), function ℛ becomes
Implementation of Algorithm 1
We now apply Algorithm 1 to determine the region in which the stable periodic orbit persists.
Step 1. We calculate the maximum value of δ such that ℛ is positive in the interior of B_{ δ }(K_{0}), in which
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
Let us denote by an equilibrium of the system, that is, f(x_{ e }(K), K) = 0. From the 7th and 6th differential equations of the model in equation (7), we have
and therefore
Substituting expression (22) in the first differential equation of (7), we have the following result
or
Equation (23) and f(x_{ e }(K), K) = 0 imply that
For the nonzero equilibrium, from equation (24) we have that
From the 3rd differential equation of (7), equation (22) and equation (24), we obtain that
From the 5th differential equation of (7) and equation (26), we have that
Substituting equations (27) and (28) in the 4th differential equation of the equation (7), we have that
Employing equation (29) in equations (21), (22), (27) and (28) we have that
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 .
To check the monotonicity of the function ℛ with respect to k_{1} over the box B_{δ*}(K_{0}), we first define a new function H as follows
in which K = [k_{1}, ⋯, k_{14}]. Then, we define the following optimization problem
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}.
The maximum value of ℛ in the box B_{δ*}(K_{0}) is computed by solving the following optimization problem
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}.
2 We consider the following optimization problem
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 perform this onedimensional parameter study, we define a path on which ℛ changes from zero to its maximum value in the box B_{δ*}(K_{0}) monotonically. Therefore, we first determine the maximum value of ℛ in the box B_{δ*}(K_{0}) by employing HGA and SQP. The maximum value of ℛ is achieved at = k_{0,1} + δ_{1}k_{0,1}, ⋯, k_{0,26} + δ_{26}k_{0,26} with
and it is 18665. We vary the value of ℛ from zero at to its maximum at by varying α from zero to one in (1  α) + α. Figure 4 shows the variation of ℛ with respect to α.
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.
References
 1.
Barkai N, Leibler S: "Robustness in simple biochemical networks". Nature. 1997, 387: 913917. 10.1038/43199
 2.
Alon U, Surette M, Barkai N, Leibler S: "Robustness in bacterial chemotaxis". Nature. 1999, 397: 168171. 10.1038/16483
 3.
von Dassow G, Meir E, Munro EM, Odell GM: "The segment polarity network is a robust development module". Nature. 2000, 406: 188192. 10.1038/35018085
 4.
Savageau MA: " Parameter sensitivity as a criterion for evaluating and comparing the performance of biochemical systems". Nature. 1971, 229: 542544. 10.1038/229542a0
 5.
Del Vecchio D: "Design and Analysis of an ActivatorRepressor Clock in E. coli". Proceedings of American Control Conference, New York. 2007
 6.
El Samad H, Del Vecchio D, Khammash M: "Repressilators and Promotilators: Loop Dynamics in Gene Regulatory Networks". Proceedings of American Control Conference, Santa Barbara. 2005
 7.
Elowitz MB, Leibler S: "A Synthetic Oscillatory Network of Transcriptional Regulators". Nature. 2000, 403: 335338. 10.1038/35002125
 8.
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.
 9.
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:20050024
 10.
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/bti348
 11.
Ma L, Iglesias PA: "Quantifying robustnes of biochemical network models". BMC Bioinformatics. 2002, 3: 3850. 10.1186/14712105338
 12.
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:20070008
 13.
Wolf J, BeckerWeiman S, Heinrich R: "Analysing the robustness of cellular rhythms". IEE Proceedings Systems Biology. 2005, 2 (1): 3541. 10.1049/sb:20045035
 14.
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.
 15.
Liu JS: "Monte Carlo strategies in scientific computing". 2001, New York: Springer
 16.
Lobo FG, Goldberg DL: "Decision making in a hybrid genetic algorithm". 1996, Technical Report No. 96009, Illinois Genetic Algorithms Laboratory
 17.
Wiggins S: "Introduction to Applied Nonlinear Dynamical Systems and Chaos". 2003, New York: Springer
 18.
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. 1964
 19.
Laub MT, Loomis WF: "A molecular network that produces spantaneous oscillations in excitable cells of Dictyostelium". Mol Biol Cell. 1998, 9: 35213532.
 20.
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)748334
 21.
Voit EO: "Computational Analysis of Biochemical System". 2000, Cambridge University Press, Cambridge, UK
 22.
Kelley CT: "Solving Nonlinear Equations with Newton's Method". Fundamentals of Algorithms, SIAM, Philadelphia. 2003
 23.
Burden RL, Faires JD: "Numerical Analysis". 2000, Thomson Brooks/cole; Florence, KY
 24.
Flecher R: "Practical methods of Optimization". 1981, 2: Constrained Optimization, Wiley, Chichester, England
 25.
Petty RH: "Neutrophil oscillations: temporal and spatiotemporal aspects of cell behaviour". Immunol Res. 2000, 23: 8594. 10.1385/IR:23:1:85.
 26.
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.1481
 27.
Michael LR, Torczon V: "Pattern Search Methods for Linearly Constrained Minimization". SIAM Journal on Optimization. 2000, 10 (3): 916941.
Acknowledgements
This work was in part supported by The Rackham School at University of Michigan.
Author information
Additional information
Authors' contributions
RG developed the computational techniques outlined in the paper, applied the techniques to the two systems described in the paper, and contributed to the writing of the paper. PI contributed to interpreting the results of the analysis on the first model, suggested the second model as one for which the proposed techniques would apply, and contributed to interpreting the results of the analysis on the second model. JS contributed to the writing and organization of the paper, and helped to enhance the numerical robustness analysis tools used in this paper. DDV developed the theoretical framework on which the proposed computational techniques are based and lead the writing and organization of the paper. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Periodic Orbit
 Hopf Bifurcation
 Oscillatory Behavior
 Sequential Quadratic Programming
 Negative Real Part