Document Detail

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:  PubMed-not-MEDLINE    
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 two-dimensional 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 isoflavone-containing soy products and endothelial function: a bayesian met...
21396465 - Responsiveness and reliability of mri in knee osteoarthritis: a meta-analysis of publis...
Publication Detail:
Type:  Journal Article     Date:  2008-11-19
Journal Detail:
Title:  EURASIP journal on bioinformatics & systems biology     Volume:  -     ISSN:  1687-4145     ISO Abbreviation:  EURASIP J Bioinform Syst Biol     Publication Date:  2009  
Date Detail:
Created Date:  2008-12-16     Completed Date:  2010-06-28     Revised Date:  2011-10-17    
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 16-18, Leipzig, Germany. nicole.radde@imise.uni-leipzig.de
Export Citation:
APA/MLA Format     Download EndNote     Download BibTex
MeSH Terms
Descriptor/Qualifier:
Comments/Corrections

From MEDLINE®/PubMed®, a database of the U.S. National Library of Medicine

Full Text
Journal Information
Journal ID (nlm-ta): EURASIP J Bioinform Syst Biol
Journal ID (publisher-id): BSB
ISSN: 1687-4145
ISSN: 1687-4153
Publisher: Hindawi Publishing Corporation
Article Information
Download PDF
Copyright © 2009 Nicole Radde.
open-access: 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: 2009E-location 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 16-18, 04107 Leipzig, Germany
Correspondence: *Nicole Radde: nicole.radde@imise.uni-leipzig.de
[other] Recommended by Matthias Steinfath

1. Introduction

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., [14] 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 switch-like behavior between two or more stable steady states [6, 810]. 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 [1315], 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 [1820]. 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. Time-scale 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 time-scale differences and time delays on a two-component 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 two-component systems, interesting recent studies indicate the impact of multiple interlocked feedback loops for the robustness of periodic behavior [3035].

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 two-component 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 two-dimensional 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 two-gene 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., [3740]), 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.


2. Stabilizing Oscillations with Time Delays
2.1. Modeling Biological Oscillators

We consider the following ODE model:

[Formula ID: eq1]
(1) 
x˙(t) = f(x(t)), x ∈ ℝ+n
with a continuously differentiable function f : ℝ + n→ℝ + n. The function f is characterized by the monotonicity of each component fi(xj) with respect to xj. This condition assigns each regulator j of i either a purely activating or a purely inhibiting function. In the first case, ∂x˙i(t)/∂xj > 0 independent of the state x of the system, the second case corresponds to ∂x˙i(t)/∂xj < 0 for all states x. Such systems can be illustrated by directed graphs G(V, E) with sign-labeled edges, often denoted interaction graphs, or, equivalently, by the signed Jacobian matrix Sf with elements sij ∈ { +, −, 0} according to the edge signs in G(V, E). Properties of this model class and the role of feedback loops are discussed in [1, 2, 4].

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 first-order 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 xs, which can be done via investigating the eigenvectors of the Jacobian matrix Jf(xs). Elements of Jf(xs) will be denoted by aij := ∂ fi(x)/∂xj|x=xs throughout the manuscript, dropping their dependence on the coordinates of xs. According to Lyapunov's indirect method [44], a hyperbolic steady-state xs is stable if all eigenvalues λ of Jf(xs) 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.

2.2. Stabilizing Oscillations with Time Delays in Two-Dimensional Systems

We consider the following delay differential equation (DDE) system:

[Formula ID: eq2]
(2) 
x˙1(t) = f1(x1(t),x2(t−τ)),  x˙2(t) = f2(x1(t),x2(t))
with monotonicity and boundedness constraints as defined above. For short, we will use the common notation xi for xi(t) and xiτ for xi(tτ) subsequently. System 2 is infinite dimensional, since the initial conditions are real-valued functions xiinit(t) : [− τ, 0]→ℝ+, which can complicate the analysis considerably. The stability of fixed points of (2) can, however, analogous to ODE systems, be determined by investigating the signs of the real parts of the roots of the characteristic equation. This characteristic equation for a fixed point xs of system (2) is derived via linearizing about xs (see also [4648]) as follows:
[Formula ID: eq3]
(3) 
χτ(λ) = det[Jf1(xs)+e−λτJf2(xs)︸=:Jfτ(xs)−λI]
with
[Formula ID: eq4]
(4) 
Jf1(xs)=(a110a21a22)  Jf2(xs) = (0a12τ00),a12τ=(∂f1(x1,x2τ)∂x2τ)|x=xs.
Calculating the determinant, χτ(λ) is given by
[Formula ID: eq6]
(5) 
χτ(λ) = (a11−λ)(a22−λ)︸=:g(λ)−e−λτa12τa21︸=:hτ(λ).
Equation (5) is a polynomial of degree two for τ = 0, which has two complex conjugate solutions λ1,2. For τ ≠ 0, it is a transcendental equation with a countable infinite number of roots. However, the number of roots in the right-half plane is known to be finite [49]. Here, we will only investigate the course of the two solutions λ1,2 in dependence of τ.

In two-component 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 two-dimensional oscillators, the activator-inhibitor oscillator (AIO) and the substrate-depletion oscillator (SDO). Both are characterized by the course of their nullclines and, related to that, the signed Jacobian matrix Sf(xs) at the unstable fixed point xs in the interior of the limit cycle, which is

[Formula ID: eq7]
(6) 
SfAIO(xs) = (+−+−) or SfSDO(xs) = (++−−)
for the AIO and the SDO, respectively. Examples for these two oscillator models can be found in [6, 26, 27]. A typical course of the nullclines for an AIO is shown in Figure 1. Note that both the AIO and the SDO do not strictly fulfill the monotonicity condition. Due to the nonlinear autocatalytic activation, the sign of the corresponding diagonal element in Jf(x) depends on the state x, as can be seen at nullcline 1 in the figure. This is, however, not a problem here, since relevant statements still hold if the state space is partitioned into regions in which the signs of the Jacobian matrix are constant [2]. Moreover, in many biological networks, a self-regulation is often not a direct interaction, but comprises intermediate components. An example is a protein that promotes expression of its own gene as a transcription factor. Here, strict monotonicity can again be reconstructed by introducing two separate variables for mRNA and protein concentration, respectively.

A stable limit cycle surrounds an unstable fixed point that is given as the intersection of the two nullclines. A fixed point xs 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 a11 > 0 at the fixed point xs. It has been shown that this condition becomes sufficient for sufficiently large time-scale 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 xs with signed Jacobian matrix of the forms SfAIO or SfSDO (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, a11 > 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 xs through a time delay)

Assume system (2) to have a stable fixed point xs forτ = 0 and signed Jacobian matrix SfAIO(xs) as given in (6). An increase in the time delay τ eventually destabilizes xs.

Theorem 1 implies that there exists a threshold time delay τth such that xs is unstable if τ > τth.

Corollary 1

If system (2) has a stable fixed point xs forτ = 0 anda11 > 0, there exists a threshold time delay τth such that xs 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 xs 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 a11 = 7.3, a22 = −27.3 and a12τa21 = −250 in (5). With these parameters, the Jacobian matrix Jfτ=0(xs) 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 xs becomes unstable through a Hopf bifurcation when the real part crosses the x-axis at τth. Further increasing τ, the imaginary parts vanish again at τ*, and both real eigenvalues eventually approach 0 and a11, respectively. In Figure 3, a12τa21 was changed to –350, and Jfτ=0(xs) 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 a11, 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 xs is exactly destabilized when they cross the imaginary axis.

The inclusion of a sufficiently large time delay τ like in (2) destabilizes xs via a Hopf bifurcation. This guarantees the existence of a stable limit cycle at least around the bifurcation value τth. Moreover, in case that xs is the only fixed point of the system, destabilizing xs 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 two-dimensional system is either a steady-state 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.

2.3. Generalizations for Higher Dimensions

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 two-dimensional 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 two-dimensional 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 xn to variable x1. The corresponding system of differential equations is given by

[Formula ID: eq8]
(7) 
x˙1= f1(x1,xnτ),x˙2= f2(x1,x2)⋮x˙n= fn(xn−1,xn)
with monotonicity and boundedness constraints as specified in Section 2. Gouzé [2] has shown that vector fields of regulatory systems with interaction graphs lacking positive circuits are injective, which implies that they have at most one fixed point. The existence of such a fixed point for the single negative loop system is shown in [53]. All together, this system has a unique fixed point, and according to [51], a stable limit cycle exists in case that this fixed point is unstable. This allows once again for a reduction of the analysis of the whole system onto the stability of its fixed point via a linearization about xs and the spectrum of Jfτ(xs).

The following theorem generalizes results in [27], where the same statement was proven for two-dimensional 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,

[Formula ID: eq9]
(8) 
dRe(λ)(τ)dτ|λ=iv > 0.
The proof can be found in Appendix C.

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 two-dimensional 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 two-dimensional systems. Unlike in two-dimensional 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.


3. On the Impact of Bifurcations for the Inverse Problem

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(ω, 𝒟):

[Formula ID: eq10]
(9) 
ω^ := arg min  ω∈ΩF(ω,𝒟).
Optimization problems of this kind are important for all fields in which differential equations are used to describe dynamic behaviors, and model parameters are to be adapted to experimental data. For nonlinear systems, these problems are known to be difficult to solve, since the surface of the objective function F has some undesirable properties. Efficient optimization algorithms are required in order to obtain reliable estimates within an acceptable time. Several approaches have been proposed (see, e.g., [39] and the subsequent discussions for an overview, or [37, 38] for a method called “multiple shooting” and applications to biological systems). Generally, these optimization problems will become even more important in systems' biology in the future.

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]:

[Formula ID: eq11]
(10) 
x˙1 = s1 − γ1x1 + kx12x12+111+x2,
[Formula ID: eq12]
(11) 
x˙2 = s2 − γx2θ+x12
with parameter vector ω = (s1, s2, γ, k, θ). The true vector is given by ω* = (0.04, 0.2, 0.035, 10, 0.1). The system has a globally stable limit cycle for these values.

We investigate two objective functions F(ω, 𝒟), first, the sum of squared errors between measurements x˜i(t) and model predictions xi(t):

[Formula ID: eq13]
(12) 
FMSE(ω,𝒟) = ∑i=1n ∑t=1T∥x˜i(t)−xi(t)∥2,
and second, the sum of differences between x˜i(t) and xi(t), weighted by the time difference Δt,
[Formula ID: eq14]
(13) 
FArea(ω,𝒟) = ∑i=1n ∑t=1TΔt⋅∥x˜i(t)−xi(t)∥.
In the limit Δt→0, (13) corresponds to the area between the two “curves” xi(t) and x˜i(t). Values for xi(t) are numerically calculated by a simple Euler discretization, and initial values are set to xi(0) = x˜i(0). By the way, numerical integration has to be performed in each step of a gradient-based optimization approach and is usually the limiting factor concerning computing time.

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.

3.1. Nonsmoothness

Figure 7 shows the bifurcation diagram of system (11) with control parameter s1. All bifurcation plots were created with the program xppaut [54]. For s1 = 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 saddle-node (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 s1 < 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 s1 = 0.02. Increasing s1, both objective functions FMSE and FArea remain almost constant. Near the bifurcation value SN, a slight increase of s1 causes an abrupt change in the qualitative dynamic behavior from convergence to oscillations (s1 = 0.03). This goes along with a jump in the objective function at the saddle-node bifurcation, which reaches zero at the true value s1* and increases smoothly thereafter.

Such jumps in the error function are generally related to bifurcations at which stable ω-limit sets disappear, typically saddle-node or subcritical Hopf bifurcations in our context. As a consequence, gradient-based methods might be much more efficient when the step size is adapted during the gradient descent to optimize F(ω, 𝒟).

3.2. Local Suboptimal Minima

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 saddle-node bifurcation (SN), where an unstable and a globally stable fixed point emerges. While the dependence of the x1-coordinate of this fixed point is only marginal (Figure 11(a)), the coordinates of x2,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 saddle-node 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.

3.3. Ruggedness Near the True Parameter Value

Figure 13 shows the bifurcation diagram with control parameter s2. 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.


4. Conclusions

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 two-dimensional systems, in particular, activator-inhibitor oscillator and substrate-depletion 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 single-loop 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 long-term 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 two-dimensional systems. Similar theorems exist for higher-dimensional systems with special interaction graphs. The single-loop 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 “mixed-circuit 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 time-consuming, 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 ill-posed 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.


Appendices
A. Proof of Theorem 1

We prove this statement by showing that one of the eigenvalues of the Jacobian matrix Jfτ(xs) approaches the positive element a11 in the limit τ, that is, lim τλ(τ) > 0 or, equivalently, lim τχτ(λ = a11) = 0. The function g(λ) in (5) is a parabola with two zeros a11 and a22, f(λ = a11) = f(λ = a22) = 0. Furthermore,

[Formula ID: eq1A.1]
(A.1) 
lim τ→∞hτ(λ) ={−∞,λ < 0,a12τa21,λ = 0,0,λ > 0.
This implies lim τhτ(λ = a11) = 0 and hence lim τχτ(λ = a11) = f(λ = a11) − lim τhτ(λ = a11) = 0, which implies further that λ = a11 > 0 is an eigenvalue of Jfτ(xs) in the limit τ. Consequently, xs cannot be stable.

B. Proof of Corollary 1

We use the following two statements for AIOs, which are derived in [27].

  1. The real part of an eigenvalue λ of the Jacobian matrix Jfτ(xs) is a continuous function of τ, Re(λ)(τ) ∈ 𝒞1(ℝ+→ℝ).
  2. 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 s1, both objective functions FMSE and FArea, 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 xs for τ = 0 whose real parts of the eigenvalues of Jf(xs) 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.

C. Proof of Theorem 2

Linearizing system (7) about xs, the characteristic equation is

[Formula ID: eq1C.1]
(C.1) 
χτ(λ) =∏i=1n(aii−λ)︸=:g(λ)−e−τλa1,nτ∏i=1nai,i−1︸=:hτ(λ)= 0.
The equation g(λ) = hτ(λ) has to be solved. For this, it is convenient to take the logarithm of both sides, which gives the condition
[Formula ID: eq1C.2]
(C.2) 
∑i=1nln (aii−λ) = −λτ + ln (a1,nτ∏i=1nai,i−1).
Calculating the derivatives of both sides with respect to λ leads to
[Formula ID: eq1C.3]
(C.3) 
d ln g(λ)dλ=−∑i=1n(aii−λ)−1,d ln hτ(λ)dλ=∂ ln g(λ)∂τdτdλ + ∂ ln hτ(λ)∂λ =−λdτdλ−τ.
Resolving for dτ/dλ yields
[Formula ID: eq1C.5]
(C.4) 
dτdλ = λ−1(∑i=1n(aii−λ)−1−τ).
Here it is worth noting that although τ is a real number, the derivative dτ/dλ is complex. It can be understood as the inverse of the derivative dλ/dτ = dRe(λ)/dτ + i(dIm(λ)/dτ), which is defined by the condition
[Formula ID: eq1C.6]
(C.5) 
(dλdτ)τ*⋅(dτdλ)λ*=λ(τ*)= 1.
We consider the real parts of both sides in (C.4):
[Formula ID: eq1C.7]
(C.6) 
Re[dτdλ] = Re[λ−1(∑i=1n(aii−λ)−τ)]
with
[Formula ID: eq1C.9]
(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[λ¯].
Inserting λ = iv and resolving for dRe[λ]/dτ show that this derivative is positive for purely imaginary eigenvalues λ:
[Formula ID: eq1C.10]
(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é J-L. 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 positive-feedback 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 enzyme-driven 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 Y-K,Cho K-H. 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. Cell-cycle control of gene expression in budding and fission yeastAnnual Review of Genetics 2005;39:69–94.
23. Chen KC,Calzone L,Csikasz-Nagy 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 time-delays on the structural stability of oscillations in a two-gene networkAdvances in Complex Systems 2008;11(3):471–483.
28. Burić N,Todorović D. Dynamics of delay-differential 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ía-Ojalvo J,Lee KH. PERIOD-TIMELESS interval timer may require an additional feedback loopPLoS Computational Biology 2007;3(8):p. e154.
32. Locke JC,Southern MM,Kozma-Bogná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 P-C,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 two-gene 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. Balsa-Canto 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 I-C,Vinga S,Vasconcelos ATR,Voit EO,Almeida JS. Parameter optimization in S-system 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 Time-Delay 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. Mallet-Paret J,Smith HL. The Poincaré-Bendixson theorem for monotone cyclic feedback systemsJournal of Dynamics and Differential Equations 1990;2(4):367–421.
52. Mallet-Paret 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 multiple-loop 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:
  • Research Article


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