A method for determining the robustness of bio-molecular oscillator models

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 Routh-Hurwitz 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][2][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 bio-molecular 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 non-oscillatory systems. Bifurcation analysis is employed to study the sensitivity of oscillatory behavior to variations of a single parameter [11][12][13][14]. Systems whose robustness to oneparameter-at-a-time 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 oneparameter-at-a-time 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 Routh-Hurwitz 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 Routh-Hurwitz stability criterion [18]. Consider a model of a biochemical network given as in which x(t) R n is a vector whose elements are concentrations of chemical species and K R p is a vector of parameters. Let K 0 R 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 S-system 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 ⊆ R p in parameter space about the nominal parameter vector K 0 = (k 0, 1 ,..., k 0, p ) and a function ℛ: R p R, called R-function, with the following properties: For all K B δ * (K 0 ), ℛ(K) > 0; There existsK on the boundary of B δ* (K 0 ) such that ℛ (K ) = 0; There exists a path Γ ending atK in D\B δ* (K 0 ) such that for all points on the path ℛ(K) < 0.
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 throughK , along which ℛ(K) is first negative, i.e., all eigenvalues of A(K) have negative real parts. Then ℛ(K) = 0 at K =K , i.e., A(K ) 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 atK , corresponding to ℛ(K ) = 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 throughK and x e (K ) is such that the terms of order higher than three are negligible in B δ* (K 0 )). According to Hopf bifurcation theorem, at K =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 R-function ℛ(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 Routh-Hurwitz 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, and T d d d d n a n n n n n : According to the Routh-Hurwitz 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 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 4. Verify that in an open neighborhood N(K ) Step 5. Verify that the value of ℛ(K) for K on the path 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 non-positive 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 ofK 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 non-significant 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 non-linear differential equations in which the state variable x = [x 1 , ..., x 7 ] represents the concentration of the seven proteins: . 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 S-system 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 ( ,..., ) x x 1 7 , in which,   12   10 4 5 7 11 13  9 3 6 8 14 12 9 3 , , 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 non-zero 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)  Table 3.
Following Remark 1, we set ℛ(K) = −d K 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 2-dimensional 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 non-positive. 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)  L . 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., d i * , which make the constraint ℛ(K) ≤ 0 active. Therefore, at the pointK defined aŝ  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.
The optimization problem is solved for j = 0, ..., 6 using the gradient descent method. We obtain the minimum values 15 , , , , and d 1 6 , 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 d j 1 , 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 -a)K 0 + aK | a > 0} ∩ B δ* + (K 0 ), we can show that it is negative except forK 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 so-called reactive oxygen species by creating superoxide O 2 − 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 1 0 1 2 1 9 2  (15) in which the state variable x = [x 1 , ..., x 16 ] represents the concentration of the 16 species. The factor r 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  (17) in which the parameter vector K R 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 0  1  1 2  13  13 14  18  0  25 26   : [ ,  , , , ,  ,  , , ,  , , ]. L k (18) 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 non-zero equilibrium renders the matrix A R 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 Routh-Hurwitz 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 v K a K d K d K L 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: min , d d i δ 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., d i * that make the constraint ℛ(K) ≥ 0 active. Therefore, at the pointK defined aŝ , 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 a K d K d K d K  Step 4. By defining a small neighborhood aboutK , that is, B .001 (K ), we obtain that d 1 11 (K) is positive in B .001 (K ) and therefore Step 4 of Algorithm 1 is completed with N(K ) := B .001 (K ) and D := N(K ) ∩ B δ* (K 0 ).
Step 5. We evaluate ℛ(K) for K on the path Γ = {(1 -aK 0 ) + aK | a > 0} ∩ B δ* + (K 0 ) and confirm that it is negative except onK 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 random-search 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 aboutK and x e (K ) 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 one-parameter-at-a-time 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 bio-molecular 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 Routh-Hurwitz 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 bio-molecular systems: the Laub    L 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) and therefore x k k k k x 7 1 13 11 14 12 = .
Substituting expression (22) in the first differential equation of (7), we have the following result 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 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   3 Determine 0 <a ≤ 1 such that ℛ(aK 0 + (1 -a)K l ) = 0 (using methods such as bisection [23]).

Appendix 5
In order to perform this one-dimensional 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 = k 0,1 + δ 1 k 0,1 , ..., k 0,26 + δ 26 k 0,26 with 0 05911 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 and it is 18665. We vary the value of ℛ from zero atK to its maximum at K by varying a from zero to one in (1 -a)K + a K . Figure 4 shows the variation of ℛ with respect to a.
In order to simulate the system over different values of ℛ, we vary the value of a 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.