Efficient Estimation of the Robustness Region of Biological Models with Oscillatory Behavior

Robustness is an essential feature of biological systems, and any mathematical model that describes such a system should reflect this feature. Especially, persistence of oscillatory behavior is an important issue. A benchmark model for this phenomenon is the Laub-Loomis model, a nonlinear model for cAMP oscillations in Dictyostelium discoideum. This model captures the most important features of biomolecular networks oscillating at constant frequencies. Nevertheless, the robustness of its oscillatory behavior is not yet fully understood. Given a system that exhibits oscillating behavior for some set of parameters, the central question of robustness is how far the parameters may be changed, such that the qualitative behavior does not change. The determination of such a “robustness region” in parameter space is an intricate task. If the number of parameters is high, it may be also time consuming. In the literature, several methods are proposed that partially tackle this problem. For example, some methods only detect particular bifurcations, or only find a relatively small box-shaped estimate for an irregularly shaped robustness region. Here, we present an approach that is much more general, and is especially designed to be efficient for systems with a large number of parameters. As an illustration, we apply the method first to a well understood low-dimensional system, the Rosenzweig-MacArthur model. This is a predator-prey model featuring satiation of the predator. It has only two parameters and its bifurcation diagram is available in the literature. We find a good agreement with the existing knowledge about this model. When we apply the new method to the high dimensional Laub-Loomis model, we obtain a much larger robustness region than reported earlier in the literature. This clearly demonstrates the power of our method. From the results, we conclude that the biological system underlying is much more robust than was realized until now.


Introduction
It is remarkable but well-known that many biological systems are robust under vastly different conditions [1,2]. Although these systems might experience strong internal or external perturbations, e.g., through environmental changes or noise, they still operate reliably. This is, for example, observed in chemotactic behavior and patterning development [2]. Robustness is an essential feature of biological systems [3,4], and any mathematical model describing their behavior should also have this property [5]. This implies the need for an efficient tool to analyze the robustness of these models.
Here we focus on the parametric robustness of biological models that show oscillatory behavior. Oscillations are ubiquitous in biology. It is found, for example, in the pulse of the heart, the circadian rythm, and the signal transduction that involves adenosine 39,59-cyclic monophospate (cAMP) in the chemotactic of Dictyostelium discoideum [6]. The robustness of a model is determined by answering the question how far the parameters of the model could be perturbed so that the qualitative behavior of the system does not change. An example of such a change is, e.g., the transition from oscillatory behavior to a steady state equilibrium. Such a drastic transition is called a Hopf bifurcation. There are many types of bifurcations possible in dynamical biological models.
Given a so-called nominal point in parameter space for which a system has a stable periodic solution, in general a region around this point exists within which the system oscillates. We call such a region a ''robustness region'' if no bifurcation of any kind occurs in its interior and if in each point of its boundary the system undergoes some bifurcation. The type of the latter bifurcations may be of any kind. An important consequence of this definition is that the period of the oscillations varies smoothly over the robustness region. If somewhere a period doubling bifurcation (also referred to as flip bifurcation) occurs, such a dramatic change in qualitative behavior indicates that the system is no longer robust. According to our definition we meet in such a point the boundary of the robustness region. (Note that in this paper, the words flip bifurcation and period doubling bifurcation are used interchangeably.) In the literature, some methods have been proposed to analyze robustness of models with oscillatory behavior. Robustness with respect to perturbations of a single or at most two parameters simultaneously can be investigated using a bifurcation analysis package such as AUTO [7]. With this package, the boundary of the robustness region can be obtained. In most cases, however, more parameters are involved and AUTO is no longer applicable. In [8], the Structured Singular Value method (SSV) from control theory [9] was used to quantify the robustness of the Laub-Loomis model [6]. This model has an oscillatory solution for a specific set of parameter values, the so-called nominal values. It was investigated how much the nominal values might be changed before a bifurcation would occur. The authors initially claimed that the allowed maximum parameter variation is 8:3%. The work was then improved by applying a hybrid optimization method which yielded a much smaller variation of 0:6% [10]. Ghaemi et al. utilized a Routh-Hurwitz stability criterion that resulted in 0:51% variation for the Laub-Loomis model [11]. The percentage values of parameter variations that are presented in these papers suggest that all parameters have the same sensitivity. However, the model might be more sensitive to some parameters than to others [12]. Furthermore, the authors studied only the Hopf bifurcation that occurs when the stable periodic behavior turns into an equilibrium.
Here we present an alternative method to analyze the parametric robustness of biological models with stable oscillatory behavior (also referred to as ''periodic solution'' or ''limit cycle''). The method aims at finding an approximation for the whole robustness region, taking into account that the sensitivity of the system might be highly parameter dependent. The consequence is that in our approach it is not useful to report the resulting estimate in terms of a percentage of the nominal value. On the contrary, the robustness region often turns out to be far from symmetric around the nominal point. Furthermore, the present approach allows for the detection of any kind of bifurcations, and is not limited to Hopf bifurcations. Another aspect refers to dimensionality. In high-dimensional systems, an important feature of any numerical method is efficiency. Many methods suffer from the so-called ''dimensional curse'', i.e. the computing time scales exponentially with dimension. For example, if we would use a Monte-Carlo approach for estimating the shape of robustness regions, we would certainly be confronted with this limiting factor. However, the method presented here has the computational advantage that it scales linearly with the number of parameters. That is the reason that it is highly efficient for systems with a high-dimensional parameter space.
The present method is based on Floquet theory and continuation of the periodic solution. Starting from the nominal parameter set, we construct an estimate for the robustness region by scanning the parameter space in orthogonal directions. If necessary, the obtained estimate is refined by shifting the nominal point to a carefully chosen new position. We do not only focus on Hopf bifurcations, but take into account all types of bifurcation that might occur to periodic solutions, including period doubling and Neimark-Sacker bifurcations. So, also bifurcations that lead to chaotic behavior may be detected. In addition, the presented method yields extra information such as the period and the amplitude of the solution for free.
To demonstrate the ideas and power of the proposed method, we apply it to ecological and biological network models that admit stable periodic solutions: the Rosenzweig-MacArthur (RMA) model and the Laub-Loomis (LL) model. The RMA model is chosen for illustrational purposes. It is well known for its rich bifurcation pattern and serves as a test case here. It is a low dimensional system for which our method is not especially designed, but it serves as a useful check of performance. It consists of three state variables with six parameters where two of them are taken free. The LL model is a high dimensional system that consists of seven state variables with fourteen parameters. Its robustness has been already investigated in [8,10,11]. As an extra check on low dimensional systems we analyze the LL model with twelve fixed and only two parameters perturbed. Our results for two dimensional systems fully agree with those obtained with existing approaches. The results for the high dimensional LL model clearly demonstrate that the present method is a real extension of the existing approaches.

Results
The stability of a periodic solution can be verified using Floquet theory (see [13] and [14]). In this theory, the Floquet multipliers, which are the eigenvalues of the so-called monodromy-matrix, are used to indicate stability. One of the Floquet multipliers is always real and equal to 1. A necessary and sufficient condition for a periodic solution to be stable is that the modulus of the other Floquet multipliers is less than 1, i.e., they lie inside the unit circle in the complex plane. If the parameters are perturbed and one of the multipliers crosses the unit circle, the solution looses its stability and a bifurcation happens. This bifurcation can be of several types as discussed in the Material and Method section.
This suggests that in order to analyze the robustness of oscillatory behavior of a model, we only need to observe its Floquet multipliers as functions of the parameters. In the Material and Method section, we describe the details to find in an efficient way an estimate for the robustness region. Starting in a so-called nominal point in parameter space for which a stable periodic solution exists, the parameter space is scanned along orthogonal directions to detect where along these lines bifurcations occur. This yields an initial estimate of the robustness region, that is gradually improved by shifting the nominal point and varying the directions.
In the next sections, we apply our method to two biological models: the low-dimensional RMA model and the high-dimensional LL model.

Application to the Rosenzweig-MacArthur Model
The Rosenzweig-MacArthur (RMA) model is an ecological model that describes the time evolution of a predator-prey system [15]. In dimensionless form, this 3-dimensional model reads as where Here, x 1 ,x 2 ,x 3 denote the prey, predator, and top predator populations, respectively, a 1 ,a 2 ,b 1 ,b 2 are the parameters in the Michaelis-Menten functions f 1 and f 2 , and k 1 ,k 2 are death rate parameters.
The dynamical behavior of this model for the fixed coefficient values has been extensively investigated in [16][17][18] as a function of k 1 and k 2 . The resulting bifurcation diagram is depicted in Figure 1A. From this figure, we see that the limit cycle behavior of the model may experience a Hopf bifurcation, a transcritical bifurcation, or may transform into a flip bifurcation. Since there are infinitely many flip bifurcations in this bifurcation diagram, it is not possible to indicate all their positions in Figure 1A. Therefore, as a warning, we shade some areas in Figure 1C to indicate that flip bifurcations may occur somewhere in these areas. Due to infinitely many flip bifurcations, we restrict ourselves to the positions of the first period doubling bifurcations, which lie on the red curved line. We apply our method to show how an estimate is obtained for the region in Figure 1A where a stable limit cycle exists. As nominal starting point we take k 0~( k 1~0 :6,k 2~0 :008). In k 0 , the solution converges to a periodic solution with period T~120:04 as shown in Figure 2A. The corresponding Floquet multipliers are m~fm 1 ,m 2 ,m 3 g~f0:9991,{4:5654e-016,{0:2319g: We notice that the largest multiplier m 1 is indeed equal to 1 within the numerical accuracy. m 2 and m 3 lie inside the unit circle, so the limit cycle in k 0 is stable. Following the method described underneath and summarized in equations (23)-(27), we construct two orthogonal directions, v 1 and v 2 , and perturb the nominal parameter set k 0 in these directions. The direction v 1 is chosen such that the Floquet multipliers will change mostly when moving along v 1 in the (k 1 ,k 2 ) plane and v 2 is orthogonal to v 1 .
Continuation is applied along perturbation direction v 1 until points B, denoted by a green star, and F, denoted by a red star, in Figure 1B are reached. Continuation is stopped at point B because the multipliers at that point are m~f1:0000,0:9991,{1:1102e-016g: So, m 2~0 :9991&1 and this indicates that the method has successfully found a fold bifurcation. Using only Floquet multipliers, one cannot discriminate between a tangent, for which the cycle collides with a saddle limit cycle, and a Hopf bifurcation, for which the limit cycle disappears into an equilibrium. However, since in both cases the boundary of the robustness region is reached, this is not a problem at all. Just for curiosity we used AUTO to confirm that it is the latter option. Continuation is stopped at point F. It does not make sense to continue beyond this point, since the value of parameter k 2 is so small there, that it is already hardly acceptable from a biological point of view. This also manifests itself in a very long period and a highly irregular shape of the limit cycle, that gives rise to a very long computational time. An example is given in Figure 2B, where we show the time behavior of the components at point F.
When the continuation procedure is applied along direction v 2 , the method hits two bifurcation points, A and D. At point A, the mutlipliers are m~f1:0006,{1:2583e-015,{1:0003g: We notice that m 3~{ 1:0003&{1, and we conclude that the method has successfully found a flip bifurcation, which is denoted by a blue-star. Since a flip indicates a possible route to chaos and it means the end of the limit cycle, as meant in the definition of robustness, this is also a boundary of robustness. On the other hand, we detect point D as a Hopf bifurcation. Thus, we obtain region ABDF as our first, crude approximation of the robustness region of the model. Note that the orthogonality of v 1 and v 2 that leads to the axes AD and BF is not directly clear from Figure 1B, because the axes have different scales.
Next, an improvement on this initial approximation is obtained by shifting the nominal point k 0 to the midpoint of the longest axis, in this case the midpoint of AD which is denoted by k Ã 0 in Figure 1C. Applying the continuation procedure to the shifted nominal point k Ã 0 along the direction v 1 , we obtain a new axis CE. Here, point C is a Hopf bifurcation point. Just as for point F, we stop continuation in E since the value of k 2 becomes too small. Together with the previous findings, we now obtain the bigger estimating region ABCDEF, as shown in Figure 1C.
During the calculations, we simultaneously obtain a lot of information on the period and the shape of the limit cycle. In fact, this information is available along all the lines through k 0 and k Ã 0 . In Figure 1D, this info is used to draw level lines for the period. It provides a nice indication how the period behaves as a function of the parameters. Since the RMA model only serves as a lowdimensional illustration of the ideas behind the proposal estimation algorithm, we will not refine the approximation further, but rather turn to a high-dimensional example.

Application to the Laub-Loomis Model
The Laub-Loomis (LL) model [6] describes the dynamical behavior of the molecular network underlying cAMP (adenosine 39,59-cyclic monophospate) oscillation observed in population of Dyctiostelium discoideum cells. The molecular network is depicted in Figure 3.
Here, after the binding of extracellular cAMP to the surface receptor CAR1, adenylate cyclase (ACA) activates internal cAMP. When internal cAMP is accumulated, it activates protein kinase PKA. In addition, ligand-bound CAR1 also activates the MAP kinase ERK2, which is then inactivated by PKA. Therefore, ERK2 no longer inhibits the cAMP phosphodiesterase REG A. A protein phosphatase activates REG A such that REG A can hydrolyze internal cAMP, hence the concentration of internal cAMP is reduced. When the internal cAMP is hydrolyzed by REG  A, PKA activity is inhibited by its regulatory subunit, so that both ACA and ERK2 activities go up.
Based on the network above, the spontaneous oscillation in cAMP is a solution of the following model At the nominal parameter set in Table 1, which is denoted by k 0 , this system has a stable periodic solution. Thus, if we choose the initial concentrations within the basin of attraction, the solution will converge to this periodic solution, as shown in Figure 4.
Restriction to a 2-dimensional cross-section of parameter space. For illustrational purposes, we first fix 12 parameters setting them at the values in Table 1 and only vary k 2 and k 14 . In this way we deal with a two dimensional cross-section in the highdimensional parameter space. The advantage is that the results can be compared to results obtained with AUTO and in [11]. AUTO yields the robustness region given in Figure 5A. This region perfectly agrees with the region reported in [11]. However, it should be noted that the method in [11] yields a very good estimate only in the two-dimensional case. For higher dimensions, their approach leads to a much more restricted estimated region. If we would apply the more-than-twodimensions approach in [11] or in [8,10] to the present twodimensional case, we would only find the small square shaped estimate indicated in Figure 6.
Applying the algorithm in (23)-(27), we obtain two directions: v 1 , which is the most sensitive direction; and v 2 , which is orthogonal to v 1 . Along these directions, we perform the continuation procedure. This leads to our first approximation of the robustness region ABDF as shown in Figure 5B.  We notice that the initial approximation is much smaller than the real robustness region found by AUTO. We improve our approximation by shifting the nominal parameter k 0 to k Ã 0 , the midpoint of AD. When the continuation procedure is applied to the new nominal parameter k Ã 0 along direction v 2 , we find the Hopf bifurcation points C and E. Together with the first approximation, we now have obtained the larger approximation region ABCDEF, as shown in Figure 5C. As extra information, we get for free the level lines for the period as indicated in Figure 5D. The approximation could be further improved by taking more perturbation directions, but this is hardly necessary to get a very good impression of the robustness region.
Application in full-dimensional parameter space. Let us now investigate the robustness region of the Laub-Loomis model in the 14-dimensional parameter space. It will be clear that in this highdimensional case the results are hard to present visually. According to algorithm (23)-(27), we find 14 orthogonal directions fv 1 ,v 2 , . . . ,v 14 g which, for convenience, are normalized to have unit length.
By applying continuation and observing the multipliers during the continuation, we easily obtain an estimate of the 14-dimensional robustness region. This estimate is shown in Figure 7A in a dedicated form. In this figure, the range of perturbations that is allowed to maintain the stability of the limit cycle is shown by horizontal lines for each perturbation direction. There are three possibilities that we stop the continuation: one of the perturbed parameters becomes close to 0 (in the LL model, all parameters should be positive), a bifurcation is detected, or the limit cycle gets an extremely long period. In the latter case, we need too much computational time to approximate the limit cycle. If one of the parameters becomes close to 0, we denote in Figure 7 the point by The continuation is stopped at c~{1:332 because then a fold bifurcation is detected, which follows from the Floquet multipliers m~f1:0002,0:9991,{1:3091e-005+1:3186e-006i,2:7715e-006, 2:2113e-016+5:501e-016ig: At c~12:6, the system still admits a stable limit cycle behavior as shown in Figure 8, but we stop the continuation because one of the perturbed parameters becomes very close to 0, see Table 2.
In the v 7 direction, the nominal parameter can be perturbed in the range k~k 0 zcv 7 , c [ ½{1:0228,4:22: ð6Þ Continuation is stopped at c~{1:0228 because the period of the limit cycle becomes extremely long and requires too much computational time. The behavior of the period along this direction is shown in Figure 9. At c~4:22, the continuation is stopped because one of the perturbed parameters becomes very close to 0. To get still a better impression of the robustness region, we shift the nominal parameter. From the result in Figure 7A, we find that the system can be mostly perturbed in the direction of v 12 . Therefore, we shift the nominal point k 0 to the midpoint of this axis, and we denote the new nominal point by k Ã 0~k 0 z5:634v 12 . When the method is applied to k Ã 0 , we obtain the results shown in Figure 7B.
Combining the information in Figures 7A and 7B, we obtain a good impression of the robustness region of the system. Contrary to the findings in [8,10,11], we conclude that the LL model has a large robustness region with a quite irregular shape.

Discussion
An important question in the modeling of biological systems is for which parameter values the model has a stable limit cycle, since this is often the parameter range in which the model describes   Table 1). In (B), the directions start in k Ã 0~k 0 z5:634v 12 . If an end point is marked with ''D'', one of the parameters has become close to zero. If an end point marked with ''*'', the period of the limit cycle becomes extremely long. If an end point does not have mark, a fold bifurcation is detected. The lengths of the horizontal lines indicate how far this direction can be followed in negative and positive directions so that a stable limit cycle is found. All directions are normalized to have unit length. A step of, e.g., length 6 in v 13 direction in (A) means that the unit vector in this direction can be made 6 times longer before a bifurcation is detected. doi:10.1371/journal.pone.0009865.g007 reality. In the literature [8,10,11], one mostly studies this topic by analyzing the eigenvalues of the Jacobian matrix of the equilibrium points of the model. For example, if some of these eigenvalues become purely imaginary, a so-called Hopf bifurcation takes place and a limit cycle comes into existence. However, analysis of eigenvalues of a Jacobian matrix is not the most appropriate way to study this problem, since these eigenvalues yield only local information. In the present paper we have presented a method to construct an estimate for the so-called robustness region in parameter space. The approach that we follow has a global, rather than a local character. Within a robustness region the system possesses a stable limit cycle and on its boundaries the system undergoes a bifurcation. A bifurcation is a dramatic change in the system dynamics indicating that the system is no longer robust if the parameters are perturbed further. For the present method, these bifurcations may be of any type and different parts of the boundary may be connected to different bifurcations.
The present method has especially been designed to scan robustness regions of systems with a high-dimensional parameter space. Its power stems from the fact that it scales linearly with the number of parameters. This implies that it is highly efficient from a numerical point of view. The present approach is based on observing the behavior of the Floquet multipliers of the periodic solution if the systems parameters are changed. In this way, one easily detects all bifurcations that may occur to the periodic solution, such as Hopf, fold, flip, and Neimark-Sacker bifurcations, which lead to disappearance or period doubling of the periodic solution.
The method has first been tested for low-dimensional systems. It is shown that for a 2-dimensional parameter space, the results are in full agreement with those obtained by the package AUTO. Thereafter the method has been applied to a high-dimensional system, the Laub-Loomis model which has 14 parameters. In this case, the method appears to be highly efficient, indeed. Contrary to the results reported in the literature [8,10,11], the method yields an estimate that is very big and irregularly shaped. The latter means that the Laub-Loomis model is much more robust with respect to changes in one parameter than in another. The present approach yields this information and is as such an extension of the methods available in literature. In the present method, a first direction is chosen such that the Floquet multipliers will change mostly if the continuation is applied along this direction. The approach finds axes that together span the estimated region.
Since all information about the limit cycle along the used axes becomes available, it requires no extra work to present, e.g., level line plots of the period of the limit cycle. Together with the general types of bifurcation that are detected, this provides a reliable and insightful impression of the dynamical behavior of a model in a wide range of values around a nominal point.

Floquet Theory and Periodic Solution
Consider an ordinary differential equation system  where x denotes the vector of state variables and k the vector of parameters. Suppose that this system has a stable periodic solution at k~k 0 with periodic solution x~x Ã and period T.
In order to investigate the stability of the solution, we linearize around the periodic orbit x Ã and obtain where J is the Jacobian matrix of (7) with respect to its state variables x. Since x Ã is T{periodic, the Jacobian matrix J is also T{periodic. According to Floquet theory (see [13] and [14]), the fundamental solution of (8), which is a matrix that is composed of n independent solutions of (8), can be written as with P(t) T{periodic and B a constant n|n matrix. Thus, Here, C~e BT is called the monodromy matrix of the system and the eigenvalues of C are called the Floquet multipliers of the system. One of them is always real and equal to 1. A necessary and sufficient condition for the periodic solution of (8) to be stable is that the other n{1 multipliers have modulus less than 1, i.e. they lie inside the unit circle. The calculation of W is explained underneath.

Calculation of Periodic Solutions
There are many methods discussed in the literature to approximate a periodic solution. To mention some of them: finite difference method, shooting method, and Poincare map method [21]. In this paper, we use the finite difference method because of its simplicity, and a short outline of the method is given below.
Consider again the ODE system (7). With the scaling with T the period, the system reads as Now, (12) has to be solved in the time interval t [ (0,1). This time interval is discretized into Nz1 points with a uniform time step h: t 1~0 ,t 2~h , . . . ,t Nz1~N h~1: The solution of (12) at time steps t~t i and t~t iz1 are related by Using the trapezoidal rule to represent the integral, we obtain where x i~x (t i ). Since the system is periodic, x(t Nz1 )~x(t 1 ), or Therefore, we have nN algebraic equations from (14) with nNz1 unknowns: Finally, since the system that we consider is autonomous, the system is invariant to a linear shift in the time origin. To remove the arbitrariness of the phase, we specify the value of one component at t~0, for example where the value g should be within the periodic solution of x 1 (t).  Thus, at time t~t 1 we have x 1~( g,x x 1 ) withx x 1 [ R n{1 . By imposing this condition, we have nN unknowns and nN algebraic equations. Its solution can be found using, e.g., Newton's scheme, provided (16) is in the orbit of x 1 (t). The details of this method can be found in [21]. So, we obtain the periodic solution in N discretized points and the value of the period T becomes known. The full periodic solution x Ã (t) can then be obtained by integrating numerically from time t~0 to t~T.
Computing Floquet Multipliers. Let us consider the principal fundamental problem, i.e. problem (8) with now d(t) taken to be a matrix with initial values where I n is the n|n identity matrix. The Floquet multipliers of the system can then be obtained by integrating the above equation for one period, that is from t~0 to t~T. Then, the Floquet multipliers, denoted by m i , i~1,2, . . . ,n, are the eigenvalues of the matrix d(T).
Note that if we employ the same numerical technique to integrate (18) and (19), both systems can be solved simultaneously. We denote the i{th column of the matrix d by d i and solve with initial conditions Since one of the multipliers should be real and equal to 1, the approximation of the periodic solution and the Floquet multipliers are carried out iteratively. If no multipliers are close to 1, we increase the number N and solve again (14) and (21) until one of the multipliers is close to 1 within a prespecified accuracy.

Continuation Method
We start at a nominal point k 0 in parameter space, where the model has a stable limit cycle, so that the Floquet multipliers lie within the unit circle (except for one). The approach outlined here is also applicable if k 0 lies on the boundary of the robustness region. The first direction v 1 , the construction of which is described below, will then point into the robustness region. It suffices to follow that direction until the boundary at the other side is met in a point k 1 , say, and to choose as new nominal point the midpoint of k 0 and k 1 . The next step is to perturb the nominal point k 0 along n orthogonal directions v 1 ,v 2 , . . . ,v m .
To construct v 1 , we introduce the function g(k)~g(k 1 ,k 2 , . . . ,k m )~max i~2,...,n which is nothing else but the largest modulus multiplier in k that is less than 1. The gradient taking e smaller and smaller until convergence is reached.
For the first direction v 1 , we now take v 1~+ g(k 0 ). For the other perturbation directions we choose vectors that are orthogonal to v 1 and to each other. They are calculated by the Gram-Schmidt method. The set of perturbation directions is thus v 1~+ g(k 0 ),v 2 , . . . ,v m f g ð26Þ Note that the choice of v 1 is unique, but the choice of v 2 , . . . ,v m is not. However, the resulting approximate for the robustness region does not much depend on this choice, unless this region is highly irregularly shaped. To check the outcome it is recommendable to apply the method with a number of different nominal points and compare the outcomes. This will give a very good impression of the situation in parameter space. The idea is now to perturb the nominal parameters k 0 along these directions, so for direction v i , we walk along the line with c both positive and negative and check for which c we approach a bifurcation. This yields the principal axes of the estimated robustness region. An improvement of this concept is obtained by repeating this procedure but with k 0 replaced by, e.g., the center of the longest axis. This leads to a refined approximation of the full robustness region. This idea is shown in Figure 1, where the initial nominal point k 0 is shifted to k ? 0 and the direction given by the line CE has been added. By another shift or by taking extra directions, this estimate can easily be improved.

Algorithm
In Figure 11 the flow chart of the algorithm is given. In this diagram we point out in a concise way that the algorithm contains the following steps: 1) Calculate the perturbation directions v j at the nominal parameter k 0 . For v 1 , take v 1~+ g(k 0 ) using (25) and construct the other perturbation directions using the Gram-Schmidt method. 2) Calculate the periodic solution and its multipliers along the lines (27) starting from k 0 . If one or more multipliers pass the unit circle, a bifurcation has been detected. 3) If refinement is required, move the nominal point to the center of the longest axis and repeat the procedure. Also, extra directions could be chosen.