The impact of time delays on the robustness of biological oscillators and the effect of bifurcations on the inverse problem.  
Jump to Full Text  
MedLine Citation:

PMID: 19079585 Owner: NLM Status: PubMednotMEDLINE 
Abstract/OtherAbstract:

Differential equation models for biological oscillators are often not robust with respect to parameter variations. They are based on chemical reaction kinetics, and solutions typically converge to a fixed point. This behavior is in contrast to real biological oscillators, which work reliably under varying conditions. Moreover, it complicates network inference from time series data. This paper investigates differential equation models for biological oscillators from two perspectives. First, we investigate the effect of time delays on the robustness of these oscillator models. In particular, we provide sufficient conditions for a time delay to cause oscillations by destabilizing a fixed point in twodimensional systems. Moreover, we show that the inclusion of a time delay also stabilizes oscillating behavior in this way in larger networks. The second part focuses on the inverse problem of estimating model parameters from time series data. Bifurcations are related to nonsmoothness and multiple local minima of the objective function. 
Authors:

Nicole Radde 
Related Documents
:

20117115  Modeling epidemics dynamics on heterogenous networks. 19658575  Random hypergraphs and their applications. 8741975  Connectionist modeling of the recovery of language functions following brain damage. 17591175  Linear matrix inequalities approach to reconstruction of biological networks. 20709515  Exposure to isoflavonecontaining soy products and endothelial function: a bayesian met... 21396465  Responsiveness and reliability of mri in knee osteoarthritis: a metaanalysis of publis... 
Publication Detail:

Type: Journal Article Date: 20081119 
Journal Detail:

Title: EURASIP journal on bioinformatics & systems biology Volume:  ISSN: 16874145 ISO Abbreviation: EURASIP J Bioinform Syst Biol Publication Date: 2009 
Date Detail:

Created Date: 20081216 Completed Date: 20100628 Revised Date: 20111017 
Medline Journal Info:

Nlm Unique ID: 101263720 Medline TA: EURASIP J Bioinform Syst Biol Country: United States 
Other Details:

Languages: eng Pagination: 327503 Citation Subset:  
Affiliation:

Institute for Medical Informatics, Statistics and Epidemiology, University of Leipzig, Härtelstrasse 1618, Leipzig, Germany. nicole.radde@imise.unileipzig.de 
Export Citation:

APA/MLA Format Download EndNote Download BibTex 
MeSH Terms  
Descriptor/Qualifier:


Comments/Corrections 
Full Text  
Journal Information Journal ID (nlmta): EURASIP J Bioinform Syst Biol Journal ID (publisherid): BSB ISSN: 16874145 ISSN: 16874153 Publisher: Hindawi Publishing Corporation 
Article Information Download PDF Copyright © 2009 Nicole Radde. openaccess: This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Received Day: 30 Month: 5 Year: 2008 Accepted Day: 14 Month: 8 Year: 2008 Print publication date: Year: 2009 Electronic publication date: Day: 19 Month: 11 Year: 2008 Volume: 2009Elocation ID: 327503 ID: 2590636 PubMed Id: 19079585 DOI: 10.1155/2009/327503 
The Impact of Time Delays on the Robustness of Biological Oscillators and the Effect of Bifurcations on the Inverse Problem  
Nicole RaddeI1*  
Institute for Medical Informatics, Statistics and Epidemiology, University of Leipzig, Härtelstraße 1618, 04107 Leipzig, Germany 

Correspondence: *Nicole Radde: nicole.radde@imise.unileipzig.de [other] Recommended by Matthias Steinfath 
The investigation of regulation mechanisms underlying various properties of cellular networks has gained much attention in recent years. Especially interesting in this setting is the relation between the topology of a regulatory network, often referred to as wiring diagram or interaction graph, and the ability of the system to exhibit certain kinds of dynamic behaviors. It is well known that feedback control mechanisms are essential for phenomena such as hysteresis, bistability, multistationarity, and periodic behavior (see, e.g., [^{1}–^{4}] or [^{5}] for a more recent review).
While these feedback mechanisms are necessary to capture such phenomena, their existence is usually by no means sufficient. A strong nonlinearity in the Goodwin oscillator model, for example, is a very restrictive requirement for oscillations [^{6}, ^{7}]. Cooperative interaction is also needed to capture switchlike behavior between two or more stable steady states [^{6}, ^{8}–^{10}]. Thus, the qualitative behavior of the system depends considerably on the exact parameter values [^{11}]. Periodic behavior, for example, can often only be observed for a small fraction in the parameter space, which is bounded by bifurcation manifolds [^{12}].
This is in contrast to real biological systems, which exhibit their function reliably under varying external conditions and internal noise [^{13}–^{15}], raising the question how this robustness is achieved. The design principles of biological networks are assumed to be a result of a long evolutionary process [^{16}, ^{17}], during which the principles are optimized for a reliable functioning. Signaling networks, for example, have to be sensitive to signals and robust against random perturbations and internal fluctuations at the same time [^{18}–^{20}]. Many cellular oscillators such as the circadian clock have to maintain a constant period and amplitude under a wide range of different external conditions [^{21}]. The regulatory network underlying the cell cycle has to be robust against perturbations, since dysfunctions may lead to programmed cell death or to phenotypes that are not able to survive for a long time, if at all [^{22}, ^{23}].
All these biological examples investigate a property of organisms which can be described by functional robustness [^{15}]. However, the exact definition of robustness varies in all these publications, which indicates that a formalization of the concept of robustness has not yet been established [^{15}, ^{24}]. Further, this goes along with the question about which mechanisms are potentially related to such a robustness.
In this paper, we focus on the robustness of biological oscillator models with respect to varying model parameters. Timescale differences, time delays, and, related to that, feedback loops comprising a large number of interactions, have already been shown to maintain periodic behavior in chemical reaction systems (see [^{1}, ^{25}] and references therein). Scheper et al. [^{25}], for example, demonstrated the importance of nonlinear regulation and time delays on a model of the circadian oscillator. Chen and Aihara [^{26}] investigated the effect of large timescale differences and time delays on a twocomponent oscillator model. Generalizations of their results can be found in [^{27}]. A stabilization of oscillations via time delays among others has also been reported in [^{21}, ^{28}, ^{29}]. While many of the earlier studies refer to twocomponent systems, interesting recent studies indicate the impact of multiple interlocked feedback loops for the robustness of periodic behavior [^{30}–^{35}].
This work focuses on oscillations induced by including a time delay into the differential equation model. This inclusion can destabilize a stable fixed point by a Hopf bifurcation. An ordinary differential equation (ODE) model describes the cell as a homogeneous chemical reaction system, assuming that the time between cause and effect of a regulatory process can be neglected. This is of course a simplification, since time delays play a role in many regulation processes. Examples are the transport of mRNA from the nucleus to the cytoplasm, diffusion processes, especially in eukaryotic cells, or the time between binding of a transcription factor to the DNA and the corresponding change in concentration of the regulated gene product. The inclusion of such time delays into the ODE models can change the dynamic behavior of the system qualitatively. Furthermore, the period of oscillations has been shown to be crucially regulated by such a delay [^{21}].
The model class considered here is characterized by monotonicity and boundedness constraints, which will be explained in detail in Section 2.1. This class is similar to systems investigated by Kaufman et al. [^{10}] and Pigolotti et al. [^{36}]. Since the proofs rely on very weak assumptions about the differential equation system, they apply to many twocomponent oscillator models which have already extensively been studied (see, e.g., [^{6}, ^{25}, ^{26}]). In this sense, the paper generalizes some of the previous publications.
Section 2.2 shows results for twodimensional systems. In particular, sufficient conditions for the destabilization of a steady state via a time delay are introduced, which imply the existence of a stable limit cycle.
Higher dimensional feedback systems are studied in Section 2.3. For a single negative feedback, system I shows that the inclusion of a time delay can destabilize a stable fixed point through a Hopf bifurcation, implying oscillating behavior. In turn, an unstable fixed point cannot become stable through a time delay.
Section 3 elucidates the problem of robustness of oscillations from a different point of view, the inference of oscillating models from time series data. We will demonstrate on a twogene network that bifurcations complicate parameter estimation considerably. They are related to nonsmooth error functions with multiple local optima. The special focus in this study is on the bifurcations relevant for chemical oscillator models. As already pointed out by several authors (see, e.g., [^{37}–^{40}]), results emphasize that advanced parameter estimation approaches for differential equations are required in this context. Finally, conclusions and ideas for future work are provided in Section 4.
We consider the following ODE model:
(1)
x˙(t) = f(x(t)), x ∈ ℝ+n 
Moreover, we assume solutions of the system to be bounded. This is a biologically plausible assumption, but excludes simple linear models. This boundedness constraint is, for example, fulfilled for all network models which describe degradation of network components as a firstorder decay process and assume bounds for the production rate [^{27}]. Solutions of these systems have the tendency to converge to steady states [^{41}, ^{42}]. Thus, more complex behavior such as oscillations is typically caused by destabilizing this steady state via Hopf bifurcations [^{43}]. Such a bifurcation requires the existence of a negative feedback loop in G(V, E) [^{2}]. Hence we focus the analysis onto the investigation of the stability of fixed points x_{s}, which can be done via investigating the eigenvectors of the Jacobian matrix J_{f}(x_{s}). Elements of J_{f}(x_{s}) will be denoted by a_{ij} := ∂ f_{i}(x)/∂x_{j}_{x=xs} throughout the manuscript, dropping their dependence on the coordinates of x_{s}. According to Lyapunov's indirect method [^{44}], a hyperbolic steadystate x_{s} is stable if all eigenvalues λ of J_{f}(x_{s}) have negative real parts [^{44}, ^{45}]. This statement also holds for differential equations including time delays [^{46}]. Further, we will refer to conditions that imply periodic behavior when a fixed point is destabilized.
We consider the following delay differential equation (DDE) system:
(2)
x˙1(t) = f1(x1(t),x2(t−τ)), x˙2(t) = f2(x1(t),x2(t)) 
(3)
χτ(λ) = det[Jf1(xs)+e−λτJf2(xs)︸=:Jfτ(xs)−λI] 
(4)
Jf1(xs)=(a110a21a22) Jf2(xs) = (0a12τ00),a12τ=(∂f1(x1,x2τ)∂x2τ)x=xs. 
(5)
χτ(λ) = (a11−λ)(a22−λ)︸=:g(λ)−e−λτa12τa21︸=:hτ(λ). 
In twocomponent ODE systems, a single negative feedback loop is not sufficient for sustained oscillations. Additionally, one of the two components must activate itself autocatalytically ([^{50}] and [^{6}, Chapter 9]). This results in two different kinds of twodimensional oscillators, the activatorinhibitor oscillator (AIO) and the substratedepletion oscillator (SDO). Both are characterized by the course of their nullclines and, related to that, the signed Jacobian matrix S_{f}(x_{s}) at the unstable fixed point x_{s} in the interior of the limit cycle, which is
(6)
SfAIO(xs) = (+−+−) or SfSDO(xs) = (++−−) 
A stable limit cycle surrounds an unstable fixed point that is given as the intersection of the two nullclines. A fixed point x_{s} can only be unstable if it is located between the minimum A and the maximum B of nullcline 1. This corresponds to a positive element a_{11} > 0 at the fixed point x_{s}. It has been shown that this condition becomes sufficient for sufficiently large timescale differences [^{26}, ^{27}]. Furthermore, it is not required any more when including time delays, which was exemplarily shown for a specific AIO model in [^{27}].
Here, we show that a stable fixed point x_{s} with signed Jacobian matrix of the forms S_{f}^{AIO} or S_{f}^{SDO} (6) can indeed always be destabilized through a time delay τ. This result is an extension of the previous results described in [^{27}], where we proved that a bifurcation via an increase in the time delay always destabilizes a fixed point. Furthermore, this destabilization is caused by a Hopf bifurcation and therefore creates a stable limit cycle. Thus, a_{11} > 0 becomes as well a sufficient condition for the existence of a stable limit cycle in AIO models, provided that the time delay is sufficiently large. This is stated in the following Theorem 1. The proof is given in Appendix A. We remark here that this theorem can analogously be proven for SDO models.
Theorem 1 (instability of x_{s} through a time delay)
Assume system (2) to have a stable fixed point x_{s} forτ = 0 and signed Jacobian matrix S_{f}^{AIO}(x_{s}) as given in (6). An increase in the time delay τ eventually destabilizes x_{s}.
Theorem 1 implies that there exists a threshold time delay τ^{th} such that x_{s} is unstable if τ > τ^{th}.
Corollary 1
If system (2) has a stable fixed point x_{s} forτ = 0 anda_{11} > 0, there exists a threshold time delay τ^{th} such that x_{s} is unstable for τ > τ^{th}. The proof is given in Appendix B
According to Theorem 1 and Corollary 1, if an AIO model has a single fixed point that lies between the minimum A and the maximum B of the first nullcline (Figure 1), this fixed point is unstable for a sufficiently large time delay τ.
Figures 2, 3, and 4 illustrate how a stable fixed point x_{s} of an AIO is destabilized by a time delay τ. According to (5), eigenvalues λ can graphically be interpreted as (probably complex) “intersections” of the two functions g(λ) and h_{τ}(λ) (Figures 2 and 3). A sketch of the real parts of the two eigenvalues λ_{1,2} as a function of τ is shown in Figure 4. For Figure 2, we have used the functions g(λ) = (x + 10)^{2} − 300 and h_{τ}(λ) = −250e^{−τλ}, that is, parameter values a_{11} = 7.3, a_{22} = −27.3 and a_{12}^{τ}a_{21} = −250 in (5). With these parameters, the Jacobian matrix J_{f}^{τ=0}(x_{s}) has two negative real eigenvalues λ_{1} and λ_{2}. Increasing τ, these eigenvalues eventually coalesce at a value τ^{A} (Figure 4(a)) and become a pair of complex conjugates, whose real parts increase with τ. The fixed point x_{s} becomes unstable through a Hopf bifurcation when the real part crosses the xaxis at τ^{th}. Further increasing τ, the imaginary parts vanish again at τ*, and both real eigenvalues eventually approach 0 and a_{11}, respectively. In Figure 3, a_{12}^{τ}a_{21} was changed to –350, and J_{f}^{τ=0}(x_{s}) has a pair of complex conjugate eigenvalues λ, λ¯ with negative real parts. The behavior of their real parts is shown in Figure 4(b) and similar to the course in Figure 4(a).
We checked this result also numerically by separating real and imaginary parts of (5) and solving for Re(λ) and Im(λ). This was done using the Newton method iteratively for a delay interval τ ∈ [0, 0.5]. Resulting real and imaginary parts are shown in Figure 5. The course is in agreement with that in Figure 4. The two eigenvalues for τ = 0 are λ1,2 =−10±50. Initial guesses λ_{init} were set using a Monte Carlo approach. Both eigenvalues present and built a pair of complex conjugates with increasing real parts. Finally, the imaginary part vanishes again and λ_{1} and λ_{2} approach 0 and a_{11}, as described before. Simulation studies (not shown here) indicate that the two eigenvalues λ_{1,2} investigated here are the rightmost ones of the spectrum, since x_{s} is exactly destabilized when they cross the imaginary axis.
The inclusion of a sufficiently large time delay τ like in (2) destabilizes x_{s} via a Hopf bifurcation. This guarantees the existence of a stable limit cycle at least around the bifurcation value τ^{th}. Moreover, in case that x_{s} is the only fixed point of the system, destabilizing x_{s} implies global convergence to a stable limit cycle from arbitrary initial conditions. This follows from the PoincaréBendixson theorem (PBT), which states that the ωlimit set of a bounded forward trajectory of a twodimensional system is either a steadystate or a limit cycle [^{45}]. In other words, a bounded solution of such a system either converges to a fixed point or to a limit cycle.
The analysis of systems with more than two variables can be more complex, since concepts of the phase plane analysis and related theorems can not always directly be transferred to the ℝ^{n}, n ≥ 3. An example is the PBT, which is often used to show the existence of a stable limit cycle in twodimensional systems. There is no general analogous theorem for higher dimensions. At least, the Hopf bifurcation theorem [^{44}] claims the existence of such a limit cycle locally about the bifurcation parameter. However, some extensions of the PBT to differential equation systems whose flow is on a twodimensional manifold [^{45}] and regulatory systems of more than two components and special graph structures have been considered [^{51}, ^{52}]. Among these systems is the single negative loop structure, which I consider in the following. It consists of a single negative feedback loop (Figure 6). The product of edge signs has to be negative. We include a time delay in the regulation of variable x_{n} to variable x_{1}. The corresponding system of differential equations is given by
(7)
x˙1= f1(x1,xnτ),x˙2= f2(x1,x2)⋮x˙n= fn(xn−1,xn) 
The following theorem generalizes results in [^{27}], where the same statement was proven for twodimensional systems.
Theorem 2
A Hopf bifurcation caused by an increase of the time delay τ in system (7) always destabilizes a stable fixed point. The pair of complex conjugate eigenvalues crosses the imaginary axis from left to right,
(8)
dRe(λ)(τ)dτλ=iv > 0. 
Theorem 2 shows that whenever an increase in the time delay τ causes a Hopf bifurcation in the single loop system, this bifurcation destabilizes a stable fixed point and creates a limit cycle.
However, sufficient conditions for the occurrence of such a bifurcation in negative loop systems similar to that in twodimensional systems remain to be investigated in this context. Also for the single loop system, the introduction of a positive autoregulation of one of the components seems to be sufficient for the existence of such a threshold value, with the same convergence argument as for twodimensional systems. Unlike in twodimensional systems, however, it is not clear for arbitrary network structures whether the destabilization of a fixed point implies the existence of a stable limit cycle not only locally in the neighborhood of a bifurcation.
In this section, the effect of bifurcations on the inverse problem to estimate parameters from time series data is investigated. We consider the following inverse problem. Given a differential equation model x˙(t) = f(x(t), ω) and time series data 𝒟 = (x˜i(t))i=1,…,n, t=1,…,T, the task is to estimate values for the parameter vector ω by minimizing an objective function F(ω, 𝒟):
(9)
ω^ := arg min ω∈ΩF(ω,𝒟). 
Here we show the impact of bifurcations on the objective function F with a special focus on bifurcations relevant for periodic behavior in regulatory systems. Results are illustrated using the AIO model described in [^{27}]:
(10)
x˙1 = s1 − γ1x1 + kx12x12+111+x2, 
(11)
x˙2 = s2 − γx2θ+x12 
We investigate two objective functions F(ω, 𝒟), first, the sum of squared errors between measurements x˜i(t) and model predictions x_{i}(t):
(12)
FMSE(ω,𝒟) = ∑i=1n ∑t=1T∥x˜i(t)−xi(t)∥2, 
(13)
FArea(ω,𝒟) = ∑i=1n ∑t=1TΔt⋅∥x˜i(t)−xi(t)∥. 
Of course, even without noise, (12) and (13) depend on the dataset 𝒟, in particular, on the initial vector x(0) and the sampling time points. Here, we show exemplary examples for fixed initial conditions x(0) = (0, 0) and simulations over the transient and two oscillation periods. Further, for simplicity reasons, we vary only one single parameter α ∈ ω at a time, the control parameter, while the rest is fixed to the true values ω*. Thus, the measurements x˜(t) are obtained via simulations using ω*, and x(t) corresponds to simulations using ω* and a different value for the control parameter α. In order to overcome the dependence of the error functions (12) and (13) on the sequence of sampling time points, we use a very small time step Δt = 0.01, which corresponds to 10.000 Euler steps for numerical integration in the simulated range.
Figure 7 shows the bifurcation diagram of system (11) with control parameter s_{1}. All bifurcation plots were created with the program xppaut [^{54}]. For s_{1} = 0, the system has two separated ωlimit sets, a stable limit cycle and a stable fixed point. Increasing τ, the stable fixed point coalesces with an unstable one in a saddlenode (SN) bifurcation, and the limit cycle becomes globally attracting. Finally, the system undergoes a supercritical Hopf bifurcation (HB), in which the second unstable fixed point becomes globally stable and the limit cycle vanishes. For a value s_{1} < SN and initial conditions x(0) = (0, 0), the system converges to the lower stable fixed point. The difference between the two curves x(t)^{true} and x(t) is large, as shown in Figure 8 and control value s_{1} = 0.02. Increasing s_{1}, both objective functions F^{MSE} and F^{Area} remain almost constant. Near the bifurcation value SN, a slight increase of s_{1} causes an abrupt change in the qualitative dynamic behavior from convergence to oscillations (s_{1} = 0.03). This goes along with a jump in the objective function at the saddlenode bifurcation, which reaches zero at the true value s_{1}* and increases smoothly thereafter.
Such jumps in the error function are generally related to bifurcations at which stable ωlimit sets disappear, typically saddlenode or subcritical Hopf bifurcations in our context. As a consequence, gradientbased methods might be much more efficient when the step size is adapted during the gradient descent to optimize F(ω, 𝒟).
Figure 10 shows the bifurcation diagram with control parameter θ. For low values of θ, the system has a stable limit cycle around an unstable fixed point. This limit cycle vanishes at a saddlenode bifurcation (SN), where an unstable and a globally stable fixed point emerges. While the dependence of the x_{1}coordinate of this fixed point is only marginal (Figure 11(a)), the coordinates of x_{2,s} increase with increasing θ (Figure 11(b)), which leads to a local suboptimal minimum in the error functions, here at a value θ = 0.86 (Figure 12). Moreover, it can be seen that the true value θ* has a relatively small basin of attraction, which is bounded by the saddlenode bifurcation (SN). On the contrary, the basin of attraction for the local minimum at θ = 0.86, which corresponds to the converging time series, is much larger.
Hence, starting a local search method with an arbitrary initial parameter vector leads in most cases to suboptimal minima which correspond to systems that converge to a stable fixed point. These local minima render global search methods such as simulated annealing or genetic algorithms necessary, which usually require long running times [^{40}]. Thus, efficient algorithms are needed in this context.
Figure 13 shows the bifurcation diagram with control parameter s_{2}. The system has a stable limit cycle bounded by two supercritical Hopf bifurcations (HBs). Within this oscillating region, period and amplitude vary considerably with the control parameter, as indicated in the simulations in Figure 14. This dependency causes multiple local minima in the error functions (Figure 15). Thus, even if the optimization process is already started within the oscillating region in the parameter space, a simple gradient search might fail to find the true parameter value but get stuck in one of the local minima.
This emphasizes again the necessity of efficient optimization algorithms for parameter estimation of differential equation models in general.
This paper investigated the robustness of sustained oscillations in regulatory systems with respect to varying model parameters. Differential equations based on chemical reaction kinetics, which are often used for this purpose, are not always robust, and oscillations only occur in a small region of the parameter space bounded by bifurcation manifolds.
In the first part of the paper, we focused on the inclusion of time delays into the differential equations. Time delays take some time between the cause and the effect of a regulation into account, and they are known to stabilize oscillations by enlarging the region in the parameter space which correspond to periodic solutions. Since the typical behavior of the class of systems considered here is convergence to a fixed point, oscillations are usually induced via destabilizing a fixed point through a Hopf bifurcation. We investigated the stability of a fixed point in dependence of the time delay. We provided sufficient conditions for a time delay to induce oscillations in twodimensional systems, in particular, activatorinhibitor oscillator and substratedepletion oscillator models, which are the typical oscillator models in two dimensions. These conditions are graphically defined in terms of the qualitative course of the nullclines, which are usually easily accessible. Specifically, if the system has a single fixed point located between the minimum and the maximum of one of the nullclines, it can always be destabilized by a sufficiently large time delay, which implies sustained oscillations. Results are based on rather general assumptions about the underlying differential equation system, which hold for many related oscillator models.
Moreover, for singleloop systems with an arbitrary number of components we showed that a Hopf bifurcation that is caused by increasing the time delay always destabilizes a stable fixed point. The real parts of the eigenvalues of the Jacobian matrix at the fixed point change signs from negative to positive. Thus, a stable fixed point can loose stability by increasing the time delay, which leads to the existence of a stable limit cycle, but an unstable fixed point cannot become stable.
Here, the analysis of the system was done by linearizing the system about a fixed point and investigating the stability of this fixed point via the spectrum of eigenvalues. This facilitates the analysis of the longterm behavior considerably. We referred to the conditions necessary for such a destabilization of a fixed point to imply sustained oscillations. The PoincaréBendixson theorem is extremely useful in this context for twodimensional systems. Similar theorems exist for higherdimensional systems with special interaction graphs. The singleloop system considered here belongs to these systems. However, for a further generalization of these results to higher dimensional systems, the following questions remain to be investigated in the future. First, can the class of systems that have a unique fixed point be further characterized? It is already known that networks lacking a positive feedback loop have at most one fixed point. Second, how can this be further generalized to networks that also contain positive feedback circuits? Networks with only positive loops cannot have stable limit cycles. Consequently, oscillations can only occur in networks that have at least one negative loop. Contrary to negative feedback control, positive loops can lead to multiple fixed points. Hence for such “mixedcircuit networks,” it is not sufficient any more to show the existence of a fixed point. It also has to be investigated whether it is unique. However, necessary conditions for multiple fixed points in positive loop systems are also known to be very restrictive, and many of these models seem to have a unique fixed point, too. Third, in which cases does a destabilization of a fixed point lead to oscillating behavior? And forth, what are sufficient conditions for the existence of a threshold time delay τ^{th}?
The second part of this work investigated the influence of bifurcations on the inverse problem to estimate model parameters from time series data. Such bifurcations are generally related to nonsmoothness and multiple local minima of the objective function to be optimized in this setting. Although these phenomena are generally not new, the focus was on the special properties occurring in the class of oscillating models considered here. Global search methods are required to find the real optimum. Together with the numerical integration, these methods are usually extremely timeconsuming, even for small systems with only a few parameters.
In a realistic setting, the problem is even worse. First, the optimization problem is of course multidimensional. All values of model parameters have in principle to be found at a time, which renders a comprehensive search and an investigation of the whole objective function difficult. Second, the data is usually noisy and sparse, leading to illposed optimization problems. Noisy datasets also mean that the measured initial condition x˜(0) might not always be the best choice for x(0). Generally, x(0) should be included into the objective function as an additional variable that has to be optimized as well, which increases the dimension of the inverse problem even further. Moreover, the dataset might also contain missing values or unobserved variables. This raises additional problems, and estimating parameters by minimizing the residual error might fail in this context anyway. In this setting, stochastic approaches might be more convenient, since they take the noise in the dataset into account. Bayesian learning approaches, for example, allow for an appropriate regularization via prior distributions over model parameters. Concluding, the development of efficient approaches for parameter estimation in differential equation models remains a challenging research field in the future.
We prove this statement by showing that one of the eigenvalues of the Jacobian matrix J_{f}^{τ}(x_{s}) approaches the positive element a_{11} in the limit τ→∞, that is, lim _{τ→∞}λ(τ) > 0 or, equivalently, lim _{τ→∞}χ_{τ}(λ = a_{11}) = 0. The function g(λ) in (5) is a parabola with two zeros a_{11} and a_{22}, f(λ = a_{11}) = f(λ = a_{22}) = 0. Furthermore,
(A.1)
lim τ→∞hτ(λ) ={−∞,λ < 0,a12τa21,λ = 0,0,λ > 0. 
We use the following two statements for AIOs, which are derived in [^{27}].
 The real part of an eigenvalue λ of the Jacobian matrix J_{f}^{τ}(x_{s}) is a continuous function of τ, Re(λ)(τ) ∈ 𝒞^{1}(ℝ_{+}→ℝ).
 Moreover, a sign change of the real part of an eigenvalue λ caused by an increase of the time delay τ is always a change from positive to negative, and includes a pair of complex conjugate eigenvalues λ_{1,2} = u ± iv with nonzero imaginary parts v ≠ 0, that is, ∂Re(λ)/∂τ_{λ=iv} > 0. Increasing s_{1}, both objective functions F^{MSE} and F^{Area}, which are shown in Figure 9, remain almost constant.
Statement (1) allows to apply the intermediate value theorem (IVT). Starting with a stable fixed point x_{s} for τ = 0 whose real parts of the eigenvalues of J_{f}(x_{s}) are negative, Re(λ)(τ) has at least one zero according to the IVT and Theorem 1. Uniqueness of this zero follows from statement 2, since a second zero would be a transition from positive to negative, that is, ∂Re(λ)/∂τ_{λ=iv} < 0, a contradiction.
Linearizing system (7) about x_{s}, the characteristic equation is
(C.1)
χτ(λ) =∏i=1n(aii−λ)︸=:g(λ)−e−τλa1,nτ∏i=1nai,i−1︸=:hτ(λ)= 0. 
(C.2)
∑i=1nln (aii−λ) = −λτ + ln (a1,nτ∏i=1nai,i−1). 
(C.3)
d ln g(λ)dλ=−∑i=1n(aii−λ)−1,d ln hτ(λ)dλ=∂ ln g(λ)∂τdτdλ + ∂ ln hτ(λ)∂λ =−λdτdλ−τ. 
(C.4)
dτdλ = λ−1(∑i=1n(aii−λ)−1−τ). 
(C.5)
(dλdτ)τ*⋅(dτdλ)λ*=λ(τ*)= 1. 
(C.6)
Re[dτdλ] = Re[λ−1(∑i=1n(aii−λ)−τ)] 
(C.7)
Re[dτdλ]=1∥dλ/dτ∥2dRe[λ]dτ,Re[λ−1(∑i=1n(aii−λ)−τ)] = Re[λ−1∑i=1n(aii−λ)−1]−Re[τλ] = 1∥λ∥2Re[λ¯∑i=1n(aii−λ)]−τ∥λ∥2Re[λ¯]. 
(C.8)
dRe[iv]dτ= ∥d(iv)/dτ∥2v2Re[−iv∑i=1n(aii−iv)−1]= ∥d(iv)/dτ∥2v2Re[∑i=1nv2aii2+v2]= ∥d(iv)/dτ∥2aii2+v2 > 0. 
References
1.  Cinquin O,Demongeot J. Roles of positive and negative feedback in biological systemsComptes Rendus Biologies 2002;325(11):1085–1095. [pmid: 12506722] 
2.  Gouzé JL. Positive and negative circuits in dynamical systemsJournal of Biological System 1998;6(21):11–15. 
3.  Snoussi EH. Necessary conditions for multistationarity and stable periodicityJournal of Biological System 1998;6(1):3–9. 
4.  Thomas R,D'Ari R. Biological Feedback. 1990Boca Raton, Fla, USA: CRC Press; 
5.  Thieffry D. Dynamical roles of biological regulatory circuitsBriefings in Bioinformatics 2007;8(4):220–225. [pmid: 17626067] 
6.  Fall CP,Marland ES,Wagner JM,Tyson JJComputational Cell Biology 2005;20New York, NY, USA: Springer; Interdisciplinary Applied Mathematics 
7.  Goodwin BC. Oscillatory behavior in enzymatic control processesAdvances in Enzyme Regulation 1965;3:425–438. [pmid: 5861813] 
8.  Angeli D,Ferrell JE Jr.,Sontag ED. Detection of multistability, bifurcations, and hysteresis in a large class of biological positivefeedback systemsProceedings of the National Academy of Sciences of the United States of America 2004;101(7):1822–1827. [pmid: 14766974] 
9.  Craciun G,Tang Y,Feinberg M. Understanding bistability in complex enzymedriven reaction networksProceedings of the National Academy of Sciences of the United States of America 2006;103(23):8697–8702. [pmid: 16735474] 
10.  Kaufman M,Soulé C,Thomas R. A new necessary condition on interaction graphs for multistationarityJournal of Theoretical Biology 2007;248(4):675–685. [pmid: 17681551] 
11.  Eißing T,Waldherr S,Allgöwer F,Scheurich P,Bullinger E. Steady state and (bi) stability evaluation of simple protease signalling networksBiosystems 2007;90(3):591–601. [pmid: 17314003] 
12.  Goldbeter A. Computational approaches to cellular rhythmsNature 2002;420(6912):238–245. [pmid: 12432409] 
13.  Barkai N,Leibler S. Circadian clocks limited by noiseNature 2000;403(6767):267–268. [pmid: 10659837] 
14.  Kitano H. Systems biology: a brief overviewScience 2002;295(5560):1662–1664. [pmid: 11872829] 
15.  Stelling J,Sauer U,Szallasi Z,Doyle FJ,Doyle J. Robustness of cellular functionsCell 2004;118(6):675–685. [pmid: 15369668] 
16.  Kollmann M,Løvdok L,Bartholomé K,Timmer J,Sourjik V. Design principles of a bacterial signalling networkNature 2005;438(7067):504–507. [pmid: 16306993] 
17.  Ma W,Lai L,Ouyang Q,Tang C. Robustness and modular design of the Drosophila segment polarity networkMolecular Systems Biology 2006;2, article 70:1–9. 
18.  Hornung G,Barkai N. Noise propagation and signaling sensitivity in biological networks: a role for positive feedbackPLoS Computational Biology 2008;4(1):p. e8. 
19.  Komarova NL,Zou X,Nie Q,Bardwell L. A theoretical framework for specificity in cell signalingMolecular Systems Biology 2005;1, article 2005.0023:1–5. 
20.  Kwon YK,Cho KH. Quantitative analysis of robustness and fragility in biological networks based on feedback dynamicsBioinformatics 2008;24(7):987–994. [pmid: 18285369] 
21.  Wang R,Zhou T,Jing Z,Chen L. Modelling periodic oscillation of biological systems with multiple timescale networksIEE Proceedings Systems Biology 2004;1(1):71–84. [pmid: 17052117] 
22.  Bähler J. Cellcycle control of gene expression in budding and fission yeastAnnual Review of Genetics 2005;39:69–94. 
23.  Chen KC,Calzone L,CsikaszNagy A,Cross FR,Novak B,Tyson JJ. Integrative analysis of cell cycle control in budding yeastMolecular Biology of the Cell 2004;15(8):3841–3862. [pmid: 15169868] 
24.  Kitano H. Towards a theory of biological robustnessMolecular Systems Biology 2007;3, article 137 
25.  Scheper T,Klinkenberg D,Pennartz C,van Pelt J. A mathematical model for the intracellular circadian rhythm generatorThe Journal of Neuroscience 1999;19(1):40–47. [pmid: 9870936] 
26.  Chen L,Aihara K. A model of periodic oscillation for genetic regulatory systemsIEEE Transactions on Circuits and Systems I 2002;49(10):1429–1436. 
27.  Radde N. The effect of time scale differences and timedelays on the structural stability of oscillations in a twogene networkAdvances in Complex Systems 2008;11(3):471–483. 
28.  Burić N,Todorović D. Dynamics of delaydifferential equations modelling immunology of tumor growthChaos, Solitons & Fractals 2002;13(4):645–655. 
29.  Schmitz S,Loeffler M,Jones JB,Lange RD,Wichmann HE. Synchrony of bone marrow proliferation and maturation as the origin of cyclic haemopoiesisCell and Tissue Kinetics 1990;23(5):425–442. [pmid: 1700930] 
30.  Clodong S,Dühring U,Kronk L,et al. Functioning and robustness of a bacterial circadian clockMolecular Systems Biology 2007;3, article 90:1–9. 
31.  Kuczenski RS,Hong KC,GarcíaOjalvo J,Lee KH. PERIODTIMELESS interval timer may require an additional feedback loopPLoS Computational Biology 2007;3(8):p. e154. 
32.  Locke JC,Southern MM,KozmaBognár L,et al. Extension of a genetic network model by iterative experimentation and mathematical analysisMolecular Systems Biology 2005;1, article 2005.0013:1–9. 
33.  Romond PC,Rustici M,Gonze D,Goldbeter A. Alternating oscillations and chaos in a model of two coupled biochemical oscillators driving successive phases of the cell cycleAnnals of the New York Academy of Sciences 1999;879:180–193. [pmid: 10415827] 
34.  Venkatesh KV,Bhartiya S,Ruhela A. Multiple feedback loops are key to a robust dynamic performance of tryptophan regulation in Escherichia coliFEBS Letters 2004;563(1–3):234–240. [pmid: 15063755] 
35.  Wagner A. Circuit topology and the evolution of robustness in twogene circadian oscillatorsProceedings of the National Academy of Sciences of the United States of America 2005;102(33):11775–11780. [pmid: 16087882] 
36.  Pigolotti S,Krishna S,Jensen MH. Oscillation patterns in negative feedback loopsProceedings of the National Academy of Sciences of the United States of America 2007;104(16):6533–6537. [pmid: 17412833] 
37.  BalsaCanto E,Peifer M,Banga JR,Timmer J,Fleck C. Hybrid optimization method with general switching strategy for parameter estimationBMC Systems Biology 2008;2, article 26:1–9. [pmid: 18171472] 
38.  Peifer M,Timmer J. Parameter estimation in ordinary differential equations for biochemical processes using the method of multiple shootingIET Systems Biology 2007;1(2):78–88. [pmid: 17441551] 
39.  Ramsay JO,Hooker G,Campbell D,Cao J. Parameter estimation for differential equations: a generalized smoothing approachJournal of the Royal Statistical Society B 2007;69(5):741–796. 
40.  Vilela M,Chou IC,Vinga S,Vasconcelos ATR,Voit EO,Almeida JS. Parameter optimization in Ssystem modelsBMC Systems Biology 2008;2, article 35:1–13. [pmid: 18171472] 
41.  Hirsch MW. Systems of differential equations that are competitive or cooperative II: convergence almost everywhereSIAM Journal on Mathematical Analysis 1985;16(3):423–439. 
42.  Smith HL. Systems of ordinary differential equations which generate an order preserving flow. A survey of resultsSIAM Review 1988;30(1):87–113. 
43.  Schmidt H,Jacobsen EW. Linear systems approach to analysis of complex dynamic behaviours in biochemical networksIEE Proceedings Systems Biology 2004;1(1):149–158. [pmid: 17052125] 
44.  Guckenheimer J,Holmes P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields 1990;42New York, NY, USA: Springer; Applied Mathematical Sciences 
45.  Perko L. Differential Equations and Dynamical Systems. 1991New York, NY, USA: Springer; Texts in Applied Mathematics 
46.  Gu K,Kharitonov VL,Chen J. Stability of TimeDelay Systems. 2003Boston, Mass, USA: Birkhäuser; 
47.  Gopalsamy K. Stability and Oscillations in Delay Differential Equations of Population Dynamics 1992;74Dordrecht, The Netherlands: Kluwer Academic Publishers; Mathematics and Its Applications 
48.  Kuang Y. Delay Differential Equations: With Applications in Population Dynamics 1993;191San Diego, Calif, USA: Academic Press; Mathematics in Science and Engineering 
49.  Jarlebring E. Critical delays and polynomial eigenvalue problems Journal of Computational and Applied Mathematics. In press 
50.  Murray JD. Mathematical Biology: An Introduction 2002;17New York, NY, USA: Springer; Interdisciplinary Applied Mathematics 
51.  MalletParet J,Smith HL. The PoincaréBendixson theorem for monotone cyclic feedback systemsJournal of Dynamics and Differential Equations 1990;2(4):367–421. 
52.  MalletParet J,Sell GR. The PoincaréBendixson theorem for monotone cyclic feedback systems with delayJournal of Differential Equations 1996;125(2):441–489. 
53.  Mees AI,Rapp PE. Periodic metabolic systems: oscillations in multipleloop negative feedback biochemical control networksJournal of Mathematical Biology 1978;5(2):99–114. [pmid: 731136] 
54.  Ermentrout B. Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students. 2002Philadelphia, Pa, USA: SIAM; 
Article Categories:

Previous Document: The chromatin remodeler SPLAYED regulates specific stress signaling pathways.
Next Document: Synthesis and antimicrobial activity of silver citrate complexes.