Modelling periodic oscillation in gene regulatory networks by cyclic ...
Bulletin of Mathematical Biology 67 (2005) 339–367 www.elsevier.com/locate/ybulm
Modelling periodic oscillation in gene regulatory networks by cyclic feedback systems Ruiqi Wanga,1, Zhujun Jingb,1, Luonan Chenc,∗ a Osaka Sangyo University, Nakagaito 3-1-1, Daito, Osaka 574-8530, Japan b Hunan Normal University, Changsha, Hunan 410081, China c Department of Electrical Engineering and Electronics, Osaka Sangyo University, Nakagaito 3-1-1, Daito,
Osaka 574-8530, Japan Received 9 March 2004; accepted 21 July 2004
Abstract In this paper, we develop a new methodology to analyze and design periodic oscillators of biological networks, in particular gene regulatory networks with multiple genes, proteins and time delays, by using negative cyclic feedback systems. We show that negative cyclic feedback networks have no stable equilibria but stable periodic orbits when certain conditions are satisfied. Specifically, we first prove the basic properties of the biological networks composed of cyclic feedback loops, and then extend our results to general cyclic feedback network with less restriction, thereby making our theoretical analysis and design of oscillators easy to implement, even for large-scale systems. Finally, we use one circadian network formed by a period protein (PER) and per mRNA, and one biologically plausible synthetic gene network, to demonstrate the theoretical results. Since there is less restriction on the network structure, the results of this paper can be expected to apply to a wide variety of areas on modelling, analyzing and designing of biological systems. © 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved.
∗ Corresponding author.
E-mail address: [email protected] (L. Chen). 1 Also at: Academy of Mathematics and System Sciences, CAS, Beijing 100080, China. 0092-8240/$30 © 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.bulm.2004.07.005
340
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
1. Introduction Periodic oscillations have been found in many different natural or physical systems, in particular biological organisms which have rhythmic phenomena at all levels with periods ranging from less than a second to years (Goldbeter, 1995; Gonze et al., 2000; Chen and Aihara, 2002a,b; Chen et al., 2004). From both theoretical and experiment viewpoints, it is a great challenging problem in biological science to model, analyze and further predict the periodic behaviors of living organisms. One of the best studied rhythmic phenomena so far is circadian oscillations, which occur with a period of 24 h and play a key physiological role in the adaptation of living organisms to their periodically varying environment. It has long been suggested that circadian rhythms are produced by limit cycle oscillators at the molecular level (Dunlap, 1999) from the regulatory feedback loops, thereby study on topics, such as, how to model, analyze, predict and even create the periodic rhythms is essential and fundamental. Up to now, several theoretical models have been successfully developed to understand circadian phenomena in Drosophila (Goldbeter, 1995; Hall, 1998; Leloup and Goldbeter, 1998, 2000; Gonze et al., 2000), Neurospora (Crosthwaite et al., 1997), cyanobacteria (Golden et al., 1997), plants (Mikkar et al., 1995) and mammals (Tei et al., 1997), which not only well agree with the experimental data but also unravel their possible mechanisms. With the rapid advances in mathematics and experiments concerning the underlying regulatory mechanisms, more detailed theoretical models and general techniques are increasingly demanded to elucidate periodic behaviors, with the consideration of time delays that are particularly important for the eukaryotes due to timeconsuming transportation or diffusion processes of molecules between the nucleus and cytoplasm in a cell (Wang et al., 2004a,b). On the other hand, in addition to the natural systems, recent progress in genetic engineering has made the design and implementation of synthetic gene networks realistic from both theoretical and experimental viewpoints (Kobayashi et al., 2002, 2003; Chen et al., 2004), in particular for simple organisms, such as Escherichia coli and yeast. Actually, from the theoretical predictions, several simple gene networks have been experimentally constructed, e.g., genetic toggle switch (Gardner et al., 2000), repressilator (Elowitz and Leibler, 2000) and other biological circuits (Becskei and Serrano, 2000). The data in these experiments well agree with the theoretical predictions, which implies that mathematical modelling is a powerful tool for designing synthetic gene regulatory networks. Such simple models clearly represent a first step towards logical cellular control by manipulating and monitoring biological processes at the DNA level, and not only can be used as building blocks to synthesize the artificial biological systems, but also have great potential for biotechnological and therapeutic applications (Hasty et al., 2001; Kobayashi et al., 2003). For synthetic switching networks, a general design procedure for the systems with positive feedback loops (Kobayashi et al., 2002, 2003) has been recently developed. Such switching networks guarantee the stable switching states without any nonequilibrium dynamics, thereby making theoretical analysis and designing tractable even for largescale systems with time delays. However, for synthetic oscillating networks, although the Repressilator and hysteresis-based oscillators were proposed (Barbai and Leibler, 2000; Chen and Aihara, 2002a,b; Wang et al., 2004a,b), there has not yet been a general network
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
341
to theoretically ensure the existence of periodic oscillations in particular for large-scale systems with the consideration of the time delays. Mathematically, there are a tremendous number of theoretical results (Goodwin, 1965; Glass and Kauffman, 1973; Hunding, 1974; Walther, 1995; Wu and Zou, 1995; Belair et al., 1996; Campbell et al., 1999; Chen and Aihara, 2002b; Chen et al., 2004) providing the sufficient conditions of limit cycles in the framework of functional differential equations (FDEs), but mainly with a few variables or with linear or certain special structures when time delays are considered. Generally, it is difficult to guarantee a system converging to a limit cycle or a sustained oscillation even for a simple-structured nonlinear system. Therefore, many important physiological factors such as time delays are simply ignored in order to reduce dimensionality and complexity of the systems. It is well known, however, that such factors may play important roles in dynamics of the biological systems. Recently, based on monotone dynamical systems (Smith, 1991, 1995), Mallet-Paret and Sell (1996a,b) introduced a discrete Lyapunov functional and successfully developed a general theory to show the existence of the omega-limit set by obtaining a Morse decomposition of the global attractor for a cyclic feedback system with time delays, which opened the door to a general inquiry into not only the topological structure but also the sufficient conditions of the existence for a specific attractor, in particular periodic attractor. In this paper, we further extend the results of Mallet-Paret and Sell (1996a,b) to more general networks and provide the sufficient conditions of limit cycles for gene regulatory systems according to both local and global analysis. Specifically, this paper aims to develop a new methodology to analyze and design biological oscillating networks with time delays, by using monotone cyclic feedback systems (Mallet-Paret and Sell, 1996a,b). One desirable property for a monotone cyclic feedback system is that there are the omega-limit sets composed of only periodic orbits and equilibria (Mallet-Paret and Sell, 1996a,b). Such a property drastically simplifies theoretical analysis and design of oscillators. We prove that negative monotone cyclic feedback systems with certain conditions have no stable equilibria but stable periodic oscillations, depending on the total time delay. Such a property is clearly ideal for designing or modelling biological oscillators. In this paper, we first give a local analysis of the oscillation network by local bifurcation theory. On the basis of the local result, we derive our main theorems, which ensure negative cyclic feedback systems to converge to periodic orbits globally. Furthermore, we extend our results to more general networks, which can have any types of interactions between nodes and are much easy to implement. Finally, to demonstrate the theoretical results, we use two examples, one from a circadian network (Goldbeter, 1995) formed by a period protein (PER) and per mRNA, and the other from a synthetic gene regulatory network, i.e., Repressilator (Elowitz and Leibler, 2000), for numerical simulations. Detailed proofs of all the theoretical results are given in Appendix A. 2. Notation and assumptions In this section, we first model a general biological network by functional differential equations, and then define the cyclic feedback network (Mallet-Paret and Sell, 1996a,b)
342
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
from several general assumptions, which are used to derive the main theoretical results in this paper. Let R+ be the set of nonnegative real numbers. Assume that a biological network is composed of n chemical components, which are proteins, mRNAs, chemical complexes, or different states of the same protein, and proteins at different locations in a cell. Then the network can be written by the following functional differential equations x(t) ˙ = f (x t ),
(1) R+n
where x(t) = (x 1 (t), . . . , x n (t)) ∈ X ⊂ is the concentrations of all components at time t ∈ R. Let C+ ≡ C([−r, 0], R+n ), where C([−r, 0], R+n ) is the space of continuous maps on [−r, 0] into R+n . x t ∈ C+ is defined by x t ≡ x(t + θ ), −r ≤ θ ≤ 0. f = ( f 1 , . . . , fn ) : C+ → R+n is the reaction rates of components, continuously differentiable, and maps a bounded subset of C+ to a bounded subset of R+n . In addition, we define N = {1, . . . , n}. Note that the reaction rates f include both synthesis and degradation rates of components. Let Xˆ ⊂ C+ be an induced space of X on [−r, 0] into R+n , i.e., φ ∈ Xˆ means φ(θ ) ∈ X for −r ≤ θ ≤ 0. A function x(t; φ) ∈ R+n is said to be a solution of Eq. (1) if it satisfies Eq. (1) for all t ≥ t0 with x(t0 + θ ; φ) = φ(θ ), −r ≤ θ ≤ 0, where φ ∈ C+ is a given initial function. To emphasize on the initial function, we define x t (φ) ≡ x(t + θ ; φ) with x t0 (φ) = x(t0 + θ ; φ) = φ(θ ), −r ≤ θ ≤ 0. For Eq. (1), orbit, equilibria, periodic orbit, omega and alpha limit sets are defined in the following ways. Definition 1 (Orbit). The orbit of Eq. (1) for the initial condition φ ∈ C+ is O+ (φ) {x t (φ) : t ≥ t0 }.
(2)
Let xˆ be the constant function equal to x for all values of its argument, i.e., x(θ ˆ )=x where −r ≤ θ ≤ 0. In other words, xˆ is the natural inclusion from x ∈ R+n to xˆ ∈ C+ by x(θ ˆ ) = x, −r ≤ θ ≤ 0. Definition 2 (Equilibria). The set of equilibria for Eq. (1) is defined by E {φ ∈ C+ : φ = xˆ for some x ∈ R+n satisfying f (x) ˆ = 0}. Definition 3 (Omega and Alpha Limit Sets). The omega limit set is defined by ω(φ) {x t (φ) : t ≥ s},
(3)
(4)
s≥0
whereas the alpha limit set is defined by α(φ) {x t (φ) : t ≤ s}.
(5)
s≤0
Definition 4 (Periodic Orbit). The orbit O+ (φ) is said to be a T-periodic orbit if x T +t (φ) = x t (φ) for all t and minimal T > 0. Next, we describe the cyclic feedback networks by mainly following the definitions and assumptions of Mallet-Paret and Sell (1996a,b), Kobayashi et al. (2003) and Wang et al. (2004a,b).
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
343
Assumption 2.1 (Monotone System). ∂ f i /∂ x j > 0, 0 (or 2, then this path is said to be a feedback loop. In addition, this feedback loop is said to be positive (or negative) if l−1 m=1 s pm+1 pm = 1 (or −1). In Fig. 2(a), there are many feedback loops, e.g., 1 → 2 → 1, 2 → 3 → 2 and 1 → 2 → 3 → 2 → 1, etc., but the largest loop, i.e., 1 → 2 → 3 → 4 → 5 → 1, is unique. Notice that Definition 8 requires l > 2. Actually, when l = 2, it is a self-feedback loop, and not included in the feedback loops in Definition 8 of the paper. Each node may have a linear or nonlinear, positive or negative self-feedback loop with no time delay in the network. All self-feedback loops are omitted in Figs. 1–3. Kobayashi et al. proved that a general biological network or Eq. (1) with only positive feedback loops has no dynamical attractors except stable equilibria (Kobayashi et al., 2002, 2003) when the system is bounded. In other words, except stable equilibria, there are neither stable periodic oscillations nor other nonequilibrium attractors for the network with only positive feedback loops. On the other hand, it is also indicated that a scalar negative delayed feedback system may have periodic solutions (Diekmann et al., 1995). According to Assumptions 2.1–2.3, Eq. (6) has only one largest feedback loop connecting all nodes, which indicates the nature of the whole system: positive or negative. Based on the monotone dynamical system theory and discrete Lyapunov functional, Mallet-Paret and Sell (1996a,b) obtained the Morse decomposition and derived the following theorem for the monotone cyclic feedback system (6) in the form of the Poincaré–Bendixson type theorem, when Assumptions 2.1–2.3 are satisfied. Let the natural phase space for Eq. (6) be C(K), where K = 0 τ21 ∪ 0 τ32 ∪ · · · ∪ 0 τn,n−1 ∪ 0 τ1,n ∪ N.
(11)
Theorem 2.1 (Poincaré–Bendixson Type Theorem). Assume that Assumptions 2.1–2.3 hold. Consider the system (6) satisfying Eq. (7) and being differentiable. Let x(t) be a solution of Eq. (6) on some time interval [t0 , ∞). Let ω(x) ⊆ C(K) denote the omegalimit set of this solution in the phase space C(K). Then either
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
347
Fig. 2. A simple case of the transformation procedure. All these three network structures have the same dynamical properties. (a) A feedback network with five nodes which satisfies Assumptions 2.1–2.4. (b) Types or signs of all interaction connected with the second node are changed. (c) On the basis of (b), types or signs of all interaction connected with the third node are changed. All self-feedback loops are omitted from the figures.
348
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Fig. 3. A scheme of a modified minimal model for circadian oscillations in PER and per mRNA with a negative feedback loop. Time delays are added into the original model to consider slow diffusion, transportation or signal transduction process of molecules between nucleus and cytosol.
1. ω(x) is a single nonconstant periodic orbit; or else 2. for each solution u(t) of Eq. (6) in ω(x), i.e., for solutions with u t ∈ ω(x) for all t ∈ R, we have α(u) ∪ ω(u) ⊆ E,
(12)
where α(u) and ω(u) denote the alpha- and omega-limit sets, respectively, of this solution, and where E ⊆ C(K) denotes the set of equilibria of Eq. (6). This theorem does not provide sufficient conditions for periodic orbits, but indicates that omega-limit sets of monotone cyclic feedback systems are composed of only periodic orbits and equilibria, which is a desirable property for modelling periodic oscillations of biological systems. In particular, the omega limit set is not empty due to the bounded assumption of Assumption 2.1. However, according to Assumption 2.3, the original networks are also restricted to special feedback structures, which may limit the potential applications. Motivated by Eq. (6) and Theorem 2.1, in this paper, we give sufficient conditions for a stable nonconstant periodic orbit and further extend our results to more general cyclic feedback networks, which have less restriction and are easy for applications. Eq. (6) or (1) with only positive feedback loops falls into the class of cooperative dynamical systems, which have been intensively investigated for both ODE and differential delay equations (Smith, 1991, 1995). Cooperative systems exhibit very regular behavior, e.g., typical solutions tend to equilibria in the case of autonomous systems, or to periodic solutions in the case of nonautonomous periodic systems. In other words, there are no stable periodic solutions for autonomous monotone dynamical systems with only positive
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
349
feedback loops (Smith, 1991, 1995; Gouze, 1998; Kobayashi et al., 2003). However, when there is a negative feedback loop in the monotone dynamical network, dynamics of Eq. (1) or (6) may become very complicated and possess a variety of nonequilibrium attractors, even in the autonomous cases. In particular, it can be shown that there is at most one equilibrium in X if all loops are nonpositive and the determinant of the Jacobian matrix of f is nonzero (Gouze, 1998). Therefore, to ensure the existence of the periodic solutions, we make the following important assumption. Assumption 2.4. The feedback loop connecting all nodes of the model for Eq. (6) is negative. Define the largest-loop to be the feedback loop connecting all nodes. Note that the largest-loop is unique for Eq. (6). When the largest-loop is negative, we say that Eq. (6) is a negative cyclic feedback network. In the paper, we prove that there exist attracting nonconstant periodic orbits for negative cyclic feedback networks, when the feedback for total one-direction interactions is sufficiently strong. In other words, we show that the omega-limit set is nonempty and is composed of only periodic orbits, which are quite different from those of cooperative systems. Fig. 2(a) is an example of an interaction graph with a negative largest-loop. The largest feedback loop, i.e., 1 → 2 → 3 → 4 → 5 → 1, is a negative feedback loop because of s21 s32 s43 s54 s15 = −1. According to Assumptions 2.2 and 2.3 and Fig. 1, clearly Theorem 2.1 for the negative cyclic feedback systems requires that 1. the largest-loop is negative and all self-feedback loops can have arbitrary signs; 2. all loops excluding the self-feedback loops and the largest-loop, are positive; 3. except the nth node, the interactions are opposite for the neighboring nodes, e.g., if e21 is +, then e32 must be −, as shown in Fig. 2(a). In this paper, besides sufficient conditions of the periodic orbits, we further extend Theorem 2.1 to general cyclic feedback systems, by eliminating 3 of the above conditions or Assumption 2.3. 3. Main results 3.1. Convergence to periodic solutions Cyclic nearest neighbor systems or cyclic feedback systems in the form of differential delay equations, i.e., Eq. (6), in which the coupling between neighbors possesses a monotonicity property, were extensively studied in Mallet-Paret and Sell (1996a,b). One important property of such systems is that the Poincaré–Bendixson type theorem holds for them. Therefore, an omega-limit set for such a system is either nonconstant periodic orbits or an equilibria if the system is bounded. When such systems are written in standard feedback forms, they are either positive (except self-feedback loops, all feedback loops including the largest-loop are positive) or negative (except self-feedback loops and negative largest-loop, all feedback loops are positive) feedback systems. In the case of positive feedback systems, such systems fall into the realm of the so-called cooperative systems.
350
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Orbits of cooperative systems or systems with all positive feedback loops have a strong tendency to converge to equilibria (Smith, 1995; Kobayashi et al., 2003). In contrast, for the negative feedback systems, attracting periodic orbits may exist, and their dynamics are quite different from cooperative systems. In this section, we derive our main theorems, and show that a cyclic feedback network can be used for modelling or designing a biological oscillator when the largest-loop is negative and the feedback for the total one-direction interactions is strong enough. Moreover, we further extend the original cyclic feedback networks to general cyclic feedback networks, which have less restriction. All detailed proofs for theoretical results are given in Appendix A. Since there is no limitation for the dimensionality, a bio-oscillator can be modelled or designed even by a large-scale system. Time delays do not change the location of an equilibrium in Eq. (6) but its stability. If all eigenvalues of the characteristic equation have negative real parts, then the equilibrium is stable and there is no oscillation locally. On the other hand, with change of a parameter, e.g., time delay τ , if any complex eigenvalue crosses the imaginary axis, then a stable equilibrium loses its stability with appearance of oscillation. Let x¯ ∈ X be an equilibrium of (6). Define
(13) for i, j ∈ N. Clearly, A(0) is the Jacobian matrix of f with respect where f i j = ∂∂xfij x=x¯ to x. Then the characteristic equation of Eq. (6) evaluated at the equilibrium x¯ is det(λI − A(λ)) = 0,
(14)
where I is the n × n identity matrix, and has the following form bn λn + bn−1 λn−1 + · · · + b0 + (−1)n+1 Be−λτ = 0, (15) n−1 where n is the number of nodes, B = f 1,n i=1 f i+1,i . b j for j = 0, . . . , n − 1 are functions of f k,k+1 f k+1,k and fi,i for 1 ≤ k ≤ n − 1 and i ∈ N except f 1,n , and bn = (−1)n . Total time delay is τ = i∈N τi+1,i . B represents the total feedback strength. Notice that b j for j = 0, . . . , n − 1 are functions of fk,k+1 f k+1,k for k ∈ N, which means that all effects of interactions between nodes k and k + 1 on b j disappear but on B exist if fk,k+1 is zero. Such a property is used in the proof of Theorem 3.2 in the appendix. From Eq. (15), clearly all time delays can be equivalently reduced to a single delay τ when stability of an equilibrium is considered, which is also indicated in Eq. (10) by the transformation of Eq. (8). Such a property drastically simplifies the analysis of Eq. (6) although there are multiple time delays. When the total time delay τ = 0, assume that the equilibrium x¯ is stable for Eq. (15). By changing the total time delay τ , we study the asymptotic stability of x¯ from the roots of the characteristic Eq. (15).
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
351
Since the solution exists and is bounded in X from the assumption, the omega limit set of Theorem 2.1 is nonempty, which means that all asymptotical solutions are periodic orbits provided that there is no stable equilibrium in X. Therefore, the basic idea in this paper is to destabilize all equilibria in X, which can be actually implemented mainly by linear analysis for most cases. We first state a theoretical result for the local analysis with respect to initial values, which guarantees a negative cyclic feedback system to converge to a local periodic orbit, by using the characteristic Eq. (15). Let mod (x, y) denote remainder after x divided by y. Theorem 3.1 (Local Convergence to Nontrivial Periodic Orbits). Assume that Assumptions 2.1–2.4 hold for Eq. (6), and the following equation has one nonzero real root v¯ at x¯ ∈ X, 2 2 i i−1 i i (−1) 2 bi v + (−1) 2 bi v − B 2 = 0, (16) i∈Ie
i∈Io
and further, 1 k−1 i−1 sgn (−1) 2 bk v¯ k (−1) 2 i bi v¯ i 2 2 B v¯ k∈Io i∈Io k i−2 − (−1) 2 bk v¯ k (−1) 2 i bi v¯ i = 0 k∈I i∈I ≥2 e
(17)
e,i
where Ie = {i ∈ N : mod (i, 2) = 0} ∪ {0}, Io = {i ∈ N : mod (i, 2) = 1} and x¯ is stable at τ = 0. Then there is a τ¯ , i (−1)n (−1) 2 bi v i 1 i∈Ie (18) τ¯ = arccos , v¯ B where the range of arccos is [0, π], and the system converges to a stable periodic orbit when τ is near τ¯ and τ > τ¯ . √ Theorem 3.1 implies that if there exists a τ¯ , at which, λ = jv¯ with j = −1 is a root of Eq. (15) and the derivative of the real part of the eigenvalue at λ = jv¯ is not zero, then the cyclic feedback system will converge to a stable periodic orbit when τ is near τ¯ and τ > τ¯ , for any initial conditions near x. ¯ This theorem shows that a biological network with a negative feedback loop can have periodic orbits, which can be proven by the Hopf bifurcation theorem for functional differential equations (Kolmanovskii and Myshkis, 1999). Moreover, we show that the Hopf bifurcation is supercritical due to the stable equilibrium x¯ before bifurcation, as indicated in Appendix A. Notice that there is no restriction on the self-feedback loop for cyclic feedback systems, i.e., for any i ∈ N, the interaction from the i th chemical component directly to itself may be either nonlinear or linear, and either positive or negative. Although the expressions of Theorem 3.1 seem to be complex, the conditions are easy to satisfy.
352
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Based on the local convergence of Theorem 3.1, we describe the global convergence conditions of existence for nontrivial periodic orbits. Theorem 3.2 (Global Convergence to Nontrivial Periodic Orbits). If Eq. (6) satisfies Assumptions 2.1–2.4 and the feedback for total one-direction interactions is sufficiently ∂f ∂ fi strong, i.e., ∂∂xf1n i ∂i+1 x i for those i with ∂ x i+1 = 0 at any equilibrium is sufficiently large, then there exists a τ¯ , i (−1)n (−1) 2 bi v i 1 i∈Ie (19) τ¯ = arccos , v¯ B where the range of arccos is [0, π], and the system converges to a stable periodic orbit for almost all initial conditions φ ∈ Xˆ , when τ > τ¯ . ∂ f 1 ∂ f i+1 ∂ fi i ∂ x i for all those i with ∂ x i+1 = 0 represents the total strength of one-direction ∂ xn interactions, i.e., the product of those interactions ∂ f i+1 /∂ x i with si,i+1 = 0 but si+1,i = 0. Theorem 3.2 is proven by showing that all equilibria are unstable due to the positive real part of an eigenvalue and the real part of that eigenvalue never returns to a negative region for any τ > τ¯ . Therefore, there is no stable equilibrium but stably periodic orbits in X due to the nonempty omega limit set of Eq. (6). Actually, we can prove that when p = 2 or 3, only b02 − B 2 < 0 is needed for existence of stable periodic orbits in Theorem 3.2, instead of sufficiently strong total one-direction interactions. Theorem 3.3 (Global Convergence to Nontrivial Periodic Orbits). Assume that Eq. (6) satisfies Assumptions 2.1–2.4. If det(A(0)) < 0 at all equilibria in X, then for almost all initial conditions φ ∈ Xˆ , x t (φ) converges to a stable periodic solution, where A(λ) is defined by (13). The conditions of Theorem 3.3 are strong, and ensure that the system converges to a stable periodic orbit, regardless of any nonnegative time delays. Theorem 3.3 is proven by showing that all equilibria in X are unstable for any nonnegative time delays. It is generally not easy to guarantee stable behaviors such as equilibria or periodic orbits even for a small network with a few components, due to the nonlinearity of the system. Theorem 3.2 means that if the feedback for total one-direction interactions is sufficiently strong, an attracting periodic orbit exists. In other words, the omega-limit set is nonempty and composed of only periodic orbits. Moreover, according to the proof of Theorem 3.2 in Appendix A, the total time delay τ affects the dynamics of Eq. (6), in contrast to the cooperative systems. On the other hand, Theorem 3.3 indicates that when the determinant of the Jacobian matrix is negative for all equilibria in X, the system converges to periodic orbits with any initial function φ ∈ Xˆ and any nonnegative time delays. As indicated in Kobayashi et al. (2003), if all feedback loops of a system are positive, then its orbits have a strong tendency to converge to equilibria. However, for negative cyclic feedback systems, when conditions of Theorem 3.2 or 3.3 hold, attracting periodic orbits exist and only periodic orbits constitute the omega-limit sets, which are quite different
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
353
from those of positive feedback networks. Therefore, if conditions of Theorem 3.2 or 3.3 are satisfied, Eq. (6) will inevitably converge to stable periodic orbits. In other words, negative cyclic feedback systems have ideal properties for oscillating networks and can be used to model and design bio-oscillators. Although Theorem 3.1 is a local convergence theorem, it becomes a global convergence result for any initial condition when there is at most one equilibrium in X, as stated in Corollary 3.1. Corollary 3.1 (Global Convergence to Nontrivial Periodic Orbits). Assume that all conditions of Theorem 3.1 hold. If det(A(0)) = 0 for all x in a convex set X, then the system converges to a stable periodic orbit for almost all initial condition φ ∈ Xˆ , when τ is near τ¯ and τ > τ¯ . By showing that there is at most one equilibrium that is unstable, we can prove the corollary, which is given in the appendix. 3.2. Generalized cyclic feedback systems Generally sign (feedback) conditions on the nonlinearities f i in Assumption 2.3 require that interactions are opposite for the neighboring nodes except the first and last nodes, which are difficult to satisfy for many biological systems. Next, we eliminate such restrictions by coordinate transformations. In other words, we extend the cyclic feedback systems to general ones with any type of interactions for the neighboring nodes. For Eq. (6), which satisfies Assumptions 2.1–2.3, choose any node i ∈ N and change types or signs of all interactions connected with the node-i , we denote the system obtained under such a transformation as y˙ (t) = g(yt )
(20)
where the transformation P for x = Py is defined by a matrix as follows: for the node i∈N 0 σ1 .. (21) P= , . 0
σn
with σi = −1 and σ j = 1 for all j but j = i . By substituting x = Py into Eq. (20), we have another FDE as follows: x(t) ˙ = f (x t ) ≡ Pg(Px t )
(22)
P−1
where P = is used. We can easily prove that Eq. (22) satisfies all Assumptions 2.1–2.3. Moreover, since P is a reversible and one-to-one map, Eq. (22) is qualitatively equivalent to Eq. (20). Therefore, we have the following theorem. Theorem 3.4 (General Cyclic Feedback Systems). For Eq. (6), assume that Assumptions 2.1–2.3, hold. The transformation of Eq. (21), which changes signs of all interactions connected with any node i (1 ≤ i ≤ n), does not change its dynamical properties. Theorem 3.4 means that Eq. (20), in which the signs of all interactions connected with any node i ∈ N are changed, is qualitatively equivalent to Eq. (22), which satisfies
354
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Assumptions 2.1–2.3. Moreover, it is also easy to show that such transformation does not change the type of any feedback loop, which implies that a negative cycle feedback network is still a negative cycle feedback network under this transformation. A simple case of the transformation procedure is shown in Fig. 2(a)–(c). The difference between the Fig. 2(a) and (b) is types or signs of all interactions connected with the second node. All cases in Fig. 2 are dynamically equivalent. An important property of this transformation is that we can get any combination of interactions under such consecutive transformations, which actually do not change type of each loop. By performing such transformations for each node, we can obtain all combinations of cyclic feedback networks with any type of interactions, which are all qualitatively equivalent. Therefore, Assumption 2.3 can be eliminated for the general cyclic feedback networks, by the following corollary. Corollary 3.2. If Assumptions 2.1 and 2.2 are satisfied, then Theorems 2.1 and 3.1–3.3 and Corollary 3.1 still hold for Eq. (6). In contrast to the special structures of the original cyclic feedback networks for Theorem 2.1, Corollary 3.2 means that general negative cyclic feedback systems with any type of interactions on each node can be used to design oscillators provided that the largestloop is negative. Fig. 2 shows different feedback network structures, which have the same dynamical properties. The results of general cyclic feedback networks actually are further applied to a hybrid monotone dynamical system with both positive and negative feedback loops (Wang et al., 2004b). 4. Numerical implementation 4.1. A minimal model for circadian rhythm with time delays In this section, we demonstrate our theoretical results by modifying a minimal model (Goldbeter, 1995) for circadian rhythms in the Drosophila, which uses a period protein (PER) and the period per gene. In particular, we add time delays into the original model (Goldbeter, 1995) to consider slow diffusion, transportation or signal transduction process of molecules between nucleus and cytosol. per mRNA, whose cytosolic concentration is denoted by M, is assumed to be synthesized in the nucleus and transferred to the cytosol with time delay τm for the translation of PER, where it is degraded. The rate of synthesis of PER is proportional to M. Three states of the protein are considered in this model, i.e., unphosphorylated (P0 ), mono-phosphorylated (P1 ) and bisphosphorylated (P2 ). In this model, we consider that PN which is a complex of P2 , behaves directly as a repressor. per mRNA affected from PN with time delay τ p , is synthesized in the nucleus and transferred to the cytosol, where it accumulates at a maximum rate vs ; and it is degraded by an enzyme with maximum rate vm and Michaelis constant K m . The rate of synthesis of the PER protein from M with time delay τm , is characterized by a first-order rate constant ks . Parameters Vi and K i (i = 1, . . . , 4) denote the maximum rate and Michaelis constant of the kinase(s) and phosphatase(s) involved in the reversible phosphorylation of P0 into P1 and P1 into P2 , respectively. The repressor PN is synthesized and degraded by the first-order rate constant k2 and k1 in cytoplasm, and
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
355
is transported into the nucleus to inhibit the mRNA in the Hill type with time delay τ p . The scheme of the model for circadian oscillations in PER and per mRNA is shown in Fig. 3. The time evolution of the five-variable model is governed by the following functional differential equations, in which all parameters and concentrations are defined with respect to the total cell volume, P1 (t) P0 (t) + ks M(t − τm ) − V1 P˙0 (t) = V2 K 2 + P1 (t) K 1 + P0 (t) P0 (t) P2 (t) P˙1 (t) = V1 + V4 − D1 (P1 (t)) K 1 + P0 (t) K 4 + P2 (t) P1 (t) (23) + k2 PN (t) − D2 (P2 (t)) P˙2 (t) = V3 K 3 + P1 (t) P˙N (t) = k1 P2 (t) − k2 PN (t) K In M(t) ˙ − vm M(t) = vs n K I + PNn (t − τ p ) K m + M(t) where D2 (P2 (t)) = V4
P2 (t) P2 (t) + k1 P2 (t) + vd K 4 + P2 (t) K d + P2 (t)
D1 (P1 (t)) = V2
P1 (t) P1 (t) + V3 . K 2 + P1 (t) K 3 + P1 (t)
and
The total quantity of PER protein, Pt is given by Pt = P0 + P1 + P2 + PN ,
(24)
where n denotes the degree of cooperativity, and K I is the threshold constant for repression. By using the transformation to a single time delay, i.e., Eqs. (8) and (23) can be equivalently changed to a system with only one time delay with identical forms as the original ones except that M(t − τm ) → M (t)
(25)
and PN (t − τ p ) is replaced by PN (t − τ ) in the last equation of Eq. (23), where total time delay is τ = τ p + τm . In the following discussion, we will use the new variable M in (25). All parameters are from Goldbeter (1995) with slight modifications, and are set as vs = 0.9125 µM h−1 , vm = 1.625 µM h−1 , K m = 0.5 µM, ks = 2.25 h−1 , vd = 2.5 µM h−1 , k1 = 6.5 h−1 , k2 = 1 h−1 , K I = 1 µM, K d = 0.1 µM, n = 4, K 1 = 0.65, K 2 = 2, K 3 = 1, K 4 = 2 µM, V2 = 3.95 µM h−1 , V3 = 12.5 µM h−1 , V4 = 6.25 µM h−1 , V1 = 10. τ (µM h−1 ) is given in each case next. Clearly Eq. (23) is a negative cyclic feedback network and satisfies Assumptions 2.1 and 2.2. Fig. 4 shows a case for a sustained oscillation generated by Eq. (23) with τ = 8, for which the period of the oscillation is around 24 h. Fig. 5 shows the trajectory toward a limit cycle as a projection onto the plane formed by the concentrations of per mRNA M and total PER protein Pt . Because perturbations do not change the period or amplitude in
356
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Fig. 4. A sustained oscillation generated from Eq. (23) based on the negative feedback regulation of per mRNA by the PER protein in Drosophila at τ = 8. The period of the oscillation is around 24 h.
Fig. 5. Evolution toward a limit cycle corresponding to a sustained oscillation shown in Fig. 4 at τ = 8. The arrow indicates the direction of movement along the orbit.
the long run, a limit cycle oscillation represents a particularly stable mode of the periodic behavior. Such stability holds with the robust nature of circadian clocks which have to maintain their amplitude and period in the changing environment. According to Eqs. (16) and (18) of Theorem 3.1, we can obtain v¯ = 0.3949 and τ¯ = 4.78, which are bifurcation values and also consistent with our numerical results.
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
357
Fig. 6. A bifurcation diagram with total time delay τ as a parameter. The solid line and dashed lines represent the equilibrium is stable or unstable respectively. The dash-dotted lines indicate the maximum and minimum values of M for the sustained oscillations. Limit cycles exist when τ > τ¯ , for which the equilibrium in this region is unstable.
Since Theorem 3.1 gives local bifurcation conditions for a periodic orbit, a sustained oscillation occurs after a bifurcation point by some control parameters. This situation is shown in Fig. 6, where the control parameter is the total time delay τ . At low values of τ , the system reaches a stable steady state. With the increase of τ , a bifurcation occurs at a critical τ¯ = 4.78. After τ > τ¯ , the steady state becomes unstable and a sustained oscillation appears. The amplitudes of the sustained oscillations are also shown in Fig. 6. When τ = τ¯ , we get a pair of imaginary roots λ = ±0.3949j for the characteristic equation. Moreover, from Fig. 6, the amplitudes increase with the time delay τ , which means that the time delay affects the amplitudes of the oscillations. When τ < τ¯ , the system reaches a stable steady state. For example, Fig. 7 shows that a trajectory asymptotically converges to an equilibrium in the (Pt , M ) plane at τ = 4.70 < τ¯ . Fig. 8 is the analysis of the effect for total time delay τ on the period T . In addition to amplitudes as shown in Fig. 6, the period of oscillation, namely T , increases with the total time delay τ in almost a linear way. Therefore, the total time delay τ can be viewed as a key parameter to control both amplitude and period of an oscillation in a biological system. 4.2. Repressilator with time delays We next examine a synthetic gene regulatory network, i.e., the Repressilator by genes cI, TetR and LacI, whose periodic oscillation was experimentally investigated in a prokaryote, i.e., E. coli in vivo (Elowitz and Leibler, 2000). For eukaryotes, due to slow processes of transportation, diffusion or signal transduction of molecules between nucleus and cytoplasm, the time delays may play an important role for the cellular dynamics. In this paper, we assume that the Repressilator is constructed in an eukaryote, such as yeast. The model can be expressed as (Elowitz and Leibler, 2000; Chen and Aihara, 2002a; Chen et al., 2004):
358
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Fig. 7. An asymptotical trajectory to a stable equilibrium in the (Pt , M ) plane at τ = 4.7.
Fig. 8. Total time delay τ and period of oscillations T .
α − α1 m i (t) 1 + p j (t − τ p j )n p˙ i (t) = βm i (t − τm i ) − βpi (t)
m˙ i (t) =
(26) (27)
where i and j have the following three pairs of values: (i = 1, j = 2), (i = 2, j = 3) and (i = 3, j = 1). n is the Hill coefficient. m i ∈ R+ are mRNAs, and pi ∈ R+ are proteins. τ p j and τm i are time delays. The regulation scheme of each gene is illustrated in Fig. 9.
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
359
Fig. 9. Repressilator with three mRNAs m i and three proteins pi , (i = 1, 2, 3), signs + and − on the edges mean s = 1 and −1, respectively. Self-feedback loops are omitted from the figure.
Let total time delay be τ=
3 i=1
τm i +
3
τ pi .
(28)
i=1
By the transformation Eq. (8), that is m 1 (t − τ + τ p2 ) → m 1 (t),
m 2 (t − τm 2 ) → m 2 (t), m 3 (t − τm 3 − τ p3 − τm 2 ) → m 3 (t), p1 (t − τ + τ p2 + τm 1 ) → p1 (t), p2 (t) → p2 (t),
(29)
p3 (t − τ p3 − τm 2 ) → p3 (t),
Eqs. (26) and (27) are equivalently converted into the following form α − α1 m i (t), m˙ i (t) = 1 + p j (t − τ p j )n p˙ i (t) = βm i (t) − βpi (t),
(30) (31)
for (i = 1, j = 2), (i = 2, j = 3) and (i = 3, j = 1), where τ p2 = τ and τ p1 = τ p3 = 0. Parameters are set as n = 2, β = 2, α1 = 1 and α = 2.1. Clearly, Eqs. (26) and (27) satisfy Assumptions 2.1 and 2.2. From Eqs. (16) and (18) of Theorem 3.1, we obtain the bifurcation values v¯ = 0.1969 and τ¯ = 11.5, which are also confirmed by numerical simulation. The bifurcation diagram and analysis of the effect for total time delay τ on the period T are shown in Figs. 10 and 11 respectively, which obviously have the similar features as the circadian model. It is easy to check that b0 = 8 and |B| = 8.5975, which means that b02 − B 2 < 0. Therefore, the conditions of Theorem 3.2 are also satisfied, which
360
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Fig. 10. Effects of τ on amplitude of m 1 . The solid line and dashed lines represent stable and unstable equilibria respectively. The dash-dotted lines indicate the maximum and minimum values of m 1 for the sustained oscillations.
Fig. 11. Total time delay τ and period of oscillations T .
implies that there is no stable equilibria but stable periodic solutions for Eqs. (26) and (27) when the feedback loop is sufficiently strong. 5. Conclusion This paper developed a new methodology to analyze and design biological oscillating networks with time delays, by using cyclic feedback systems. One desirable property for a monotone cyclic feedback system is that there exists a nonempty omega-limit set composed
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
361
of only periodic orbits and equilibria. Such a property drastically simplifies theoretical analysis and design of oscillators. We first provided sufficient conditions for existence of limit cycles from both local and global analysis, and then further extended the original cyclic feedback systems into general cyclic feedback networks, which allow any types of interactions on each node. Such an extension relaxes the restrictions of cyclic feedback systems, and enables our model to apply to a wide variety of systems for modelling and designing bio-oscillators. Finally, two examples are used for numerical simulation, one from a circadian network (Goldbeter, 1995) formed by a period protein (PER) and per mRNA, and the other from a synthetic gene regulatory network, i.e., Repressilator (Elowitz and Leibler, 2000). Since time delays may play an important role particularly in eukaryotes, we added time delays into the models to study their effects on the cellular dynamics. We showed that a periodic oscillation is generally produced by the negative feedback in this model. Since the time delays may significantly influence the dynamics of the overall system, we also analyzed the effects of the total delay on the stability of the equilibria and periodic oscillations. Simulation results indicated that the total time delay can slowly change the amplitudes and vary periods in almost a linear manner. Our simulation results imply that total time delay can be used as an important control parameter both quantitatively and qualitatively to change dynamics of the systems, in particular for the amplitudes and periods of periodic oscillations. Although we have mainly examined effects of time delays on the cellular dynamics, there are also other important facts, which may play crucial roles in biological processes and should be further investigated in future works from both theoretical and experimental viewpoints, such as stochastic noise (Elowitz et al., 2002; Swain et al., 2002; Zhou et al., submitted for publication) and discrete nature of biological models. Acknowledgements The work is partly supported by Osaka Sangyo University, and National Natural Science Foundation of China (10371037). Appendix A. Proofs of theorems and corollary In this appendix, we give detailed proofs of all theorems and corollary in this paper. A.1. Proof of Theorem 3.1 With the Assumptions 2.1–2.4, Eq. (1) can be written as Eq. (6), which is a monotone negative cyclic feedback system with delays. The roots of the transcendental Eq. (15) for λ determine the stability of equilibria. If real parts of all roots are negative for an equilibrium, then it is stable. On the other hand, if there exists a root with a positive real part, the equilibrium becomes unstable. √ In other words, a Hopf bifurcation occurs when λ = jv is a root of Eq. (15), where j = −1 and v
362
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
is a nonzero real number, that is i (−1) 2 bi v i = (−1)n cos(vτ )B, i∈Ie
(−1)
i−1 2
(A.1)
bi v i = (−1)n+1 sin(vτ )B,
i∈Io
where Ie = {i ∈ N : mod (i, 2) = 0} ∪ {0} and Io = {i ∈ N : mod (i, 2) = 1}. For any τ , by eliminating terms of sin and cos, we have
2 i 2
(−1) bi v
i
+
i∈Ie
2 (−1)
i−1 2
bi v
− B 2 = 0.
i
(A.2)
i∈Io
Hence, if Eq. (A.2) has a nonzero real root v, ¯ then Eq. (15) has an imaginary root λ = jv. ¯ Notice that B < 0 due to the negative cyclic feedback loop. Therefore, we can obtain the critical values for v¯ and τ¯ , which has the following form τ¯ =
1 arccos v¯
(−1)n
i
(−1) 2 bi v¯ i
i∈Ie
B
,
(A.3)
where the range of arccos is [0, π]. On the other hand, from Eq. (15), we have
dλ dτ
−1
= (−1)n+1
h (λ) λτ τ e − , Bλ λ
(A.4)
where h (λ) = (−1)n nλn−1 + (n − 1)bn−1 λn−2 + · · · + b1 .
(A.5)
Note that i−1 i−2 1 i i (−1) 2 i bi v¯ + j (−1) 2 i bi v¯ h (λ)|λ=jv¯ = v¯ i∈I i∈I ,i≥2
o
1 (h 1 + j h 2 ). v¯
e
(A.6)
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
363
λ) dλ −1 Since sgn{ d(Re dτ } = sgn Re{( dτ ) }, by substituting Eq. (A.1) into the real part of Eq. (A.4) at λ = jv, ¯ we obtain ∂(Re λ) sgn ∂τ λ=jv¯ n+1 h (λ) λτ e = sgn Re (−1) Bλ λ=jv¯ ! h sin( vτ ¯ ) + h ¯ ) 2 cos(vτ n+1 1 = sgn (−1) (A.7) B v¯ 2 k−1 i−1 1 k 2 bk v = sgn (−1) ¯ (−1) 2 i bi v¯ i B 2v¯ 2 k∈I i∈Io o k i−2 k i . − (−1) 2 bk v¯ (−1) 2 i bi v¯ k∈Ie
i∈Ie ,i≥2
Thus, when the right part of Eq. (A.7) is nonzero, we have
∂(Re λ) ∂τ λ=jv¯
> 0, which implies
that real parts of the complex eigenvalues move to the right half plane for increasing τ with λ on the imaginary axis. Then the stable equilibrium x¯ loses its stability, according to Hopf bifurcation theorem for delay equations. Therefore, Eq. (6) has a stable periodic solution when τ is near τ¯ and τ > τ¯ . In other words, the Hopf bifurcation is supercritical. In fact, by assumption, x¯ is a stable equilibrium at τ = 0, which means when 0 < τ < τ¯ , λ) < 0 means there all real parts of eigenvalues are negative. On the other hand, ∂(Re ∂τ λ=jv¯
are at least a pair of complex eigenvalues with positive real parts when 0 < τ < τ¯ , which contradicts our assumption. This completes the proof. A.2. Proof of Theorem 3.2 Assume that feedback for total one-direction interactions is sufficiently strong, i.e., ∂ fi+1 ∂ fi i ∂ x i for those i with ∂ x i+1 = 0 at any equilibrium is sufficiently large. Then, according to the definition of B from Eq. (15), |B| is also sufficiently large, which implies that ∂ f1 ∂ xn
b02 − B 2 < 0.
(A.8)
In other words, at any equilibrium, the even-order Eq. (A.2) at least has a nonzero real root v, ¯ which means that characteristic Eq. (15) has an imaginary root λ = jv. ¯ Note that B < 0 due to the assumption of the negative cyclic feedback network. On the other hand, if B > 0, Eq. (6) is actually a network with only positive feedback loops, which has no stable periodic orbit in X according to Kobayashi et al. (2003) and Smith (1995). Moreover, using Eqs. (A.1) and (A.2), we have the following relationship for sufficiently large |B|: 1
v¯ ∼ |B| n ,
(A.9)
364
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
and for k = 0, 1, 2, . . .. τ¯k =
1 2kπ + arccos v¯
(−1)n
i
(−1) 2 bi v¯ i
i∈Ie
B
,
where the range of arccos is [0, π]. From Eq. (A.4), we have −1 dλ h (λ) λτ τ e − , = (−1)n+1 dτ Bλ λ and d(Re λ) sgn dτ
! = sgn Re
dλ dτ
−1
(A.10)
(A.11)
.
Therefore, we obtain n v¯ 2n−2 d(Re λ) = sgn sgn (1 + O( )) , dτ λ=jv¯ B2
(A.12)
(A.13)
where the expression O( ) = O( n1 ((b12 −2b0b2 )v¯ 2−2n +c2n−4 v¯ 4−2n +· · ·+c4 v¯ −4 +c2 v¯ −2 )) represents a higher order of , which is small when B is sufficiently large, where ck (k = 2, 4, . . . , 2n − 4) are functions of bi , (i = 0, 1, . . . , n − 1). When the feedback for total one-direction interactions is sufficiently strong, |B| is d(Re λ) > 0 for any number of nodes. Therefore, sufficiently large. Thus we have dτ λ=jv¯
the real part of λ moves to the right half plane for τ > τ¯ when λ is on the imaginary axis, i.e., the stable equilibrium loses the stability. However, by the assumption, the solution of Eq. (1) or (6) exists and is bounded in X, which implies that there exists a nonempty omega limit set composed of only periodic orbits and equilibria according to the Poincaré–Bendixson theorem for cyclically coupled differential delay equations (Mallet-Paret and Sell, 1996a,b) or Theorem 2.1. Since such arguments hold for all equilibria, all equilibria lose their stability and the omega limit sets are composed of only stable periodic orbits, which proves this theorem. A.3. Proof of Theorem 3.3 Consider the characteristic Eq. (15) evaluated at an equilibrium x¯ det(λI − A(λ)) = 0,
(A.14)
where I is the n × n identity matrix, and has the following form h(λ) ≡ bn λn + bn−1 λn−1 + · · · + b1 λ + b0 + (−1)n+1 Be−λτ = 0,
(A.15)
where n is the number of nodes and bn = (−1)n . Assume A(0) < 0. We have h(0) = b0 − B = A(0) < 0 and h(+∞) = +∞ for an even number n; h(0) = b0 + B = −A(0) > 0 and h(+∞) = −∞ for an odd number n. By continuality of h(λ), h(λ) = 0 has at least
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
365
one positive real root for both of the above situations. Therefore, x¯ is unstable. If the above conditions hold at all equilibria, then all equilibria are unstable. By assumption, the solution of Eq. (6) exists and is bounded in X. Therefore, according to the Poincaré–Bendixson theorem for systems of cyclically coupled differential delay equations (Mallet-Paret and Sell, 1996a,b) or Theorem 2.1, there is a nonempty omega limit set composed of only stable periodic orbits, which proves this theorem. A.4. Proof of Theorem 3.4 For Eq. (6), which satisfies Assumptions 2.1–2.3, choose any node i ∈ N and change types or signs of all interactions connected with the node-i . We denote the system obtained under such a transformation as y˙ (t) = g(yt )
(A.16)
where the transformation P is defined by a matrix as follows: for the node i ∈ N 0 σ1 .. P= , . 0
(A.17)
σn
with σi = −1 and σ j = 1 for all j but j = i . By substituting x = Py into Eq. (A.16), we have another FDE as follows: x(t) ˙ = f (x t ) = Pg(Px t )
(A.18)
where P = P−1 is used. According to Eq. (A.17), the mapping from g to f is one to one and reversible. Therefore, f is topologically equivalent to g and Eq. (A.18) has the same dynamical properties as Eq. (A.16). Moreover, because of ∂ f i /∂ x i−1 = σi σi−1 ∂gi /∂yi−1 with σi = −1 and σi−1 = 1, ∂ f i /∂ x i−1 and ∂gi /∂yi−1 have opposite signs. In the same way, it is easy to check that ∂ f i+1 /∂ x i , ∂ f i−1 /∂ x i and ∂ f i /∂ x i+1 also change their signs. However, all other feedbacks, do not change their signs. Clearly Eq. (A.18) satisfies all our Assumptions 2.1 and 2.2. This completes the proof. A.5. Proof of Corollary 3.1 To prove the corollary, we only need to show that f is injective on the convex set X or there is at most one equilibrium, when A(0) is not singular for all x ∈ X according to Gouze (1998). Assume that x 1 and x 2 are two different points in the convex set X such that f (x 1 ) = f (x 2 ). Then, " 1 D f ((1 − α)x 1 + αx 2 )(x 2 − x 1 )dα 0 = f (x 1 ) − f (x 2 ) = 0
"
1
=
A(0)dα(x 2 − x 1 )
0
where D f = A(0) is the Jacobian matrix of f with respect to x.
(A.19)
366
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Clearly, if A(0) is not singular for all x ∈ X, then the determinant det(A(0)) is nonzero and has a constant sign due to continuity of the system. Since (1 − α)x 1 + αx 2 is also #1 in the convex set X for 0 ≤ α ≤ 1, det( 0 A(0)dα) is nonzero or the integrated matrix #1 0 A(0)dα is nonsingular, which implies x 1 = x 2 from Eq. (A.19). This contradicts the assumption, and proves the corollary. References Barbai, N., Leibler, S., 2000. Biological rhythms: circadian clocks limited by noise. Nature 403, 267. Becskei, A., Serrano, L., 2000. Engineering stability in gene networks by autoregulation. Nature 405, 590–593. Belair, J., Campbell, S.A., Driessche, P., 1996. Frustration, stability, and delay-induced oscillations in a neural network model. SIAM J. Appl. Math. 56, 245–255. Campbell, S.A., Ruan, S., Wei, J., 1999. Qualitative analysis of a neural network model with multiple time delays. Int. J. Bifurcation Chaos 9, 1585–1595. Chen, L., Aihara, K., 2002a. Stability of genetic regulatory networks with time delay. IEEE Trans. Circuits Syst. I. 49, 602–608. Chen, L., Aihara, K., 2002b. A model of periodic oscillation for genetic regulatory systems. IEEE Trans. Circuits Syst. I. 49, 1429–1436. Chen, L., Wang, R., Kobayashi, T., Aihara, K., 2004. Dynamics of gene regulatory networks with cell division cycles. Phys. Rev. E 70, 011909, 1–13. Crosthwaite, S.K., Dunlap, J.C., Loros, J.J., 1997. Neurospora wc-1 and wc-2: transcription, photoresponses, and the origin of circadian rhythmicity. Science 276, 763–769. Diekmann, O., van Gils, S.A., Verduyn Lunel, S.M., Walther, H.O., 1995. Delay Equations Functional-, Complex-, and Nonlinear Analysis. Springer-Verlag. Dunlap, J.C., 1999. Molecular bases for circadian clocks. Cell 96, 271–290. Elowitz, M.B., Leibler, S., 2000. A synthetic oscillatory network of transcriptional regulators. Nature 403, 335–338. Elowitz, M.B., Levine, A.J., Siggia, E.D., Swain, P.S., 2002. Stochastic gene expression in a single cell. Science 297, 1183–1186. Gardner, T.S., Cantor, C.R., Collins, J.J., 2000. Construction of a genetic toggle switch in Escherichia Coli. Nature 403, 339–342. Glass, L., Kauffman, S.A., 1973. The logical analysis of continuous, non-linear biochemical control networks. J. Theor. Biol. 39, 103–129. Goldbeter, A., 1995. A model for circadian oscillations in the Drosophila period protein (PER). Proc. R. Soc. Lond. B 261, 319–324. Golden, S.S., Ishiura, M., Johnson, C.H., Kondo, T., 1997. Cyanobacterial circadian rhythms. Ann. Rev. Plant Physiol. Plant Mol. Biol. 48, 327–354. Gonze, D., Leloup, J.-C., Goldbeter, A., 2000. Theoretical models for circadian rhythms in Neurospora and Drosophila. C. R. Acad. Sci. Paris, III 323, 57–67. Goodwin, B.C., 1965. Oscillatory behavior in enzymatic control process. Adv. Enzyme Regua. 3, 425–438. Gouze, J.-L., 1998. Positive and negative circuits in dynamical systems. J. Biol. Syst. 6, 11–15. Hall, J.C., 1998. Genetics of biological rhythms in Drosophila. Adv. Genet. 38, 135–184. Hasty, J., Isaacs, F., Dolnik, M., McMillen, D., Colins, J.J., 2001. Designer gene networks: towards fundamental cellular control. Chaos 11, 207–220. Hunding, A., 1974. Limit cycles in enzyme systems with nonlinear negative feedback. Biophys. Struct. Mech. 1, 47–54. Kobayashi, T., Chen, L., Aihara, K., 2002. Design of genetic switches with only positive feedback loops. In: Proceedings of 2002 IEEE Computer Society Bioinformatics Conference. pp. 151–162. Kobayashi, T., Chen, L., Aihara, K., 2003. Modelling genetic switches with positive feedback loops. J. Theor. Biol. 221, 379–399. Kolmanovskii, V., Myshkis, A., 1999. Introduction to the Theory and Applications of Functional Differential Equations. Kluwer Academic Publishers.
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
367
Leloup, J.-C., Goldbeter, A., 1998. A model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins. J. Biol. Rhythms 13, 70–87. Leloup, J.-C., Goldbeter, A., 2000. Modelling the molecular regulatory mechanism of circadian rhythms in Drosophila. BioEssays 22, 84–93. Mallet-Paret, J., Sell, G.R., 1996a. Systems of differential delay equations: floquet multipliers and discrete Lyapunov Functions. J. Differ. Eqns. 125, 385–440. Mallet-Paret, J., Sell, G.R., 1996b. The Poincaré–Bendixson theorem for monotone cyclic feedback systems with delay. J. Differ. Eqns. 125, 441–489. Mikkar, A.J., Carre, I.A., Strayer, C.A., Chua, N.H., Kay, S.A., 1995. Circadian clock mutants in Arabidopsis identified by iuciferase imaging. Science 267, 1161–1163. Smith, H., 1991. Periodic tridiagonal competitive and cooperative systems of differential equations. SIAM J. Math. Anal. 22, 1102–1109. Smith, H., 1995. Monotone Dynamical Systems, vol. 41. American Mathematical Society, Providence, RI. Swain, P.S., Elowitz, M.B., Siggia, E.D., 2002. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc. Natl. Acad. Sci. 99, 12795–12800. Tei, H., Okamura, H., Shigeyoshi, Y., Fukuhara, C., Ozawa, R., Hirose, M., Sakaki, Y., 1997. Circadian oscillation of a mammalian homologue of the Drosophila period gene. Nature 389, 512–516. Walther, H.O., 1995. The 2-dimensional attractor of x(t) = −µx(t) + f (x(t − 1)). Memoirs of the Amer. Math. Soc., 544. Amer. Math. Soc, Providence, RI. Wang, R., Jing, Z., Chen, L., 2004a. Periodic oscillators in genetic networks with negative feedback loops. WSEAS Trans. Math. 3, 150–156. Wang, R., Zhou, T., Jing, Z., Chen, L., 2004b. Modelling periodic oscillation of biological systems with multiple time scale networks. Syst. Biol. 1, 71–84. Wu, J., Zou, X., 1995. Patterns of sustained oscillations in neural networks with delayed interactions. Appl. Math. Comput. 73, 55–75. Zhou, T., Chen, L., Aihara, K., 2004. Intercellular communications induced by random fluctuations (submitted for publication).
Modelling periodic oscillation in gene regulatory networks by cyclic feedback systems Ruiqi Wanga,1, Zhujun Jingb,1, Luonan Chenc,∗ a Osaka Sangyo University, Nakagaito 3-1-1, Daito, Osaka 574-8530, Japan b Hunan Normal University, Changsha, Hunan 410081, China c Department of Electrical Engineering and Electronics, Osaka Sangyo University, Nakagaito 3-1-1, Daito,
Osaka 574-8530, Japan Received 9 March 2004; accepted 21 July 2004
Abstract In this paper, we develop a new methodology to analyze and design periodic oscillators of biological networks, in particular gene regulatory networks with multiple genes, proteins and time delays, by using negative cyclic feedback systems. We show that negative cyclic feedback networks have no stable equilibria but stable periodic orbits when certain conditions are satisfied. Specifically, we first prove the basic properties of the biological networks composed of cyclic feedback loops, and then extend our results to general cyclic feedback network with less restriction, thereby making our theoretical analysis and design of oscillators easy to implement, even for large-scale systems. Finally, we use one circadian network formed by a period protein (PER) and per mRNA, and one biologically plausible synthetic gene network, to demonstrate the theoretical results. Since there is less restriction on the network structure, the results of this paper can be expected to apply to a wide variety of areas on modelling, analyzing and designing of biological systems. © 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved.
∗ Corresponding author.
E-mail address: [email protected] (L. Chen). 1 Also at: Academy of Mathematics and System Sciences, CAS, Beijing 100080, China. 0092-8240/$30 © 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.bulm.2004.07.005
340
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
1. Introduction Periodic oscillations have been found in many different natural or physical systems, in particular biological organisms which have rhythmic phenomena at all levels with periods ranging from less than a second to years (Goldbeter, 1995; Gonze et al., 2000; Chen and Aihara, 2002a,b; Chen et al., 2004). From both theoretical and experiment viewpoints, it is a great challenging problem in biological science to model, analyze and further predict the periodic behaviors of living organisms. One of the best studied rhythmic phenomena so far is circadian oscillations, which occur with a period of 24 h and play a key physiological role in the adaptation of living organisms to their periodically varying environment. It has long been suggested that circadian rhythms are produced by limit cycle oscillators at the molecular level (Dunlap, 1999) from the regulatory feedback loops, thereby study on topics, such as, how to model, analyze, predict and even create the periodic rhythms is essential and fundamental. Up to now, several theoretical models have been successfully developed to understand circadian phenomena in Drosophila (Goldbeter, 1995; Hall, 1998; Leloup and Goldbeter, 1998, 2000; Gonze et al., 2000), Neurospora (Crosthwaite et al., 1997), cyanobacteria (Golden et al., 1997), plants (Mikkar et al., 1995) and mammals (Tei et al., 1997), which not only well agree with the experimental data but also unravel their possible mechanisms. With the rapid advances in mathematics and experiments concerning the underlying regulatory mechanisms, more detailed theoretical models and general techniques are increasingly demanded to elucidate periodic behaviors, with the consideration of time delays that are particularly important for the eukaryotes due to timeconsuming transportation or diffusion processes of molecules between the nucleus and cytoplasm in a cell (Wang et al., 2004a,b). On the other hand, in addition to the natural systems, recent progress in genetic engineering has made the design and implementation of synthetic gene networks realistic from both theoretical and experimental viewpoints (Kobayashi et al., 2002, 2003; Chen et al., 2004), in particular for simple organisms, such as Escherichia coli and yeast. Actually, from the theoretical predictions, several simple gene networks have been experimentally constructed, e.g., genetic toggle switch (Gardner et al., 2000), repressilator (Elowitz and Leibler, 2000) and other biological circuits (Becskei and Serrano, 2000). The data in these experiments well agree with the theoretical predictions, which implies that mathematical modelling is a powerful tool for designing synthetic gene regulatory networks. Such simple models clearly represent a first step towards logical cellular control by manipulating and monitoring biological processes at the DNA level, and not only can be used as building blocks to synthesize the artificial biological systems, but also have great potential for biotechnological and therapeutic applications (Hasty et al., 2001; Kobayashi et al., 2003). For synthetic switching networks, a general design procedure for the systems with positive feedback loops (Kobayashi et al., 2002, 2003) has been recently developed. Such switching networks guarantee the stable switching states without any nonequilibrium dynamics, thereby making theoretical analysis and designing tractable even for largescale systems with time delays. However, for synthetic oscillating networks, although the Repressilator and hysteresis-based oscillators were proposed (Barbai and Leibler, 2000; Chen and Aihara, 2002a,b; Wang et al., 2004a,b), there has not yet been a general network
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
341
to theoretically ensure the existence of periodic oscillations in particular for large-scale systems with the consideration of the time delays. Mathematically, there are a tremendous number of theoretical results (Goodwin, 1965; Glass and Kauffman, 1973; Hunding, 1974; Walther, 1995; Wu and Zou, 1995; Belair et al., 1996; Campbell et al., 1999; Chen and Aihara, 2002b; Chen et al., 2004) providing the sufficient conditions of limit cycles in the framework of functional differential equations (FDEs), but mainly with a few variables or with linear or certain special structures when time delays are considered. Generally, it is difficult to guarantee a system converging to a limit cycle or a sustained oscillation even for a simple-structured nonlinear system. Therefore, many important physiological factors such as time delays are simply ignored in order to reduce dimensionality and complexity of the systems. It is well known, however, that such factors may play important roles in dynamics of the biological systems. Recently, based on monotone dynamical systems (Smith, 1991, 1995), Mallet-Paret and Sell (1996a,b) introduced a discrete Lyapunov functional and successfully developed a general theory to show the existence of the omega-limit set by obtaining a Morse decomposition of the global attractor for a cyclic feedback system with time delays, which opened the door to a general inquiry into not only the topological structure but also the sufficient conditions of the existence for a specific attractor, in particular periodic attractor. In this paper, we further extend the results of Mallet-Paret and Sell (1996a,b) to more general networks and provide the sufficient conditions of limit cycles for gene regulatory systems according to both local and global analysis. Specifically, this paper aims to develop a new methodology to analyze and design biological oscillating networks with time delays, by using monotone cyclic feedback systems (Mallet-Paret and Sell, 1996a,b). One desirable property for a monotone cyclic feedback system is that there are the omega-limit sets composed of only periodic orbits and equilibria (Mallet-Paret and Sell, 1996a,b). Such a property drastically simplifies theoretical analysis and design of oscillators. We prove that negative monotone cyclic feedback systems with certain conditions have no stable equilibria but stable periodic oscillations, depending on the total time delay. Such a property is clearly ideal for designing or modelling biological oscillators. In this paper, we first give a local analysis of the oscillation network by local bifurcation theory. On the basis of the local result, we derive our main theorems, which ensure negative cyclic feedback systems to converge to periodic orbits globally. Furthermore, we extend our results to more general networks, which can have any types of interactions between nodes and are much easy to implement. Finally, to demonstrate the theoretical results, we use two examples, one from a circadian network (Goldbeter, 1995) formed by a period protein (PER) and per mRNA, and the other from a synthetic gene regulatory network, i.e., Repressilator (Elowitz and Leibler, 2000), for numerical simulations. Detailed proofs of all the theoretical results are given in Appendix A. 2. Notation and assumptions In this section, we first model a general biological network by functional differential equations, and then define the cyclic feedback network (Mallet-Paret and Sell, 1996a,b)
342
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
from several general assumptions, which are used to derive the main theoretical results in this paper. Let R+ be the set of nonnegative real numbers. Assume that a biological network is composed of n chemical components, which are proteins, mRNAs, chemical complexes, or different states of the same protein, and proteins at different locations in a cell. Then the network can be written by the following functional differential equations x(t) ˙ = f (x t ),
(1) R+n
where x(t) = (x 1 (t), . . . , x n (t)) ∈ X ⊂ is the concentrations of all components at time t ∈ R. Let C+ ≡ C([−r, 0], R+n ), where C([−r, 0], R+n ) is the space of continuous maps on [−r, 0] into R+n . x t ∈ C+ is defined by x t ≡ x(t + θ ), −r ≤ θ ≤ 0. f = ( f 1 , . . . , fn ) : C+ → R+n is the reaction rates of components, continuously differentiable, and maps a bounded subset of C+ to a bounded subset of R+n . In addition, we define N = {1, . . . , n}. Note that the reaction rates f include both synthesis and degradation rates of components. Let Xˆ ⊂ C+ be an induced space of X on [−r, 0] into R+n , i.e., φ ∈ Xˆ means φ(θ ) ∈ X for −r ≤ θ ≤ 0. A function x(t; φ) ∈ R+n is said to be a solution of Eq. (1) if it satisfies Eq. (1) for all t ≥ t0 with x(t0 + θ ; φ) = φ(θ ), −r ≤ θ ≤ 0, where φ ∈ C+ is a given initial function. To emphasize on the initial function, we define x t (φ) ≡ x(t + θ ; φ) with x t0 (φ) = x(t0 + θ ; φ) = φ(θ ), −r ≤ θ ≤ 0. For Eq. (1), orbit, equilibria, periodic orbit, omega and alpha limit sets are defined in the following ways. Definition 1 (Orbit). The orbit of Eq. (1) for the initial condition φ ∈ C+ is O+ (φ) {x t (φ) : t ≥ t0 }.
(2)
Let xˆ be the constant function equal to x for all values of its argument, i.e., x(θ ˆ )=x where −r ≤ θ ≤ 0. In other words, xˆ is the natural inclusion from x ∈ R+n to xˆ ∈ C+ by x(θ ˆ ) = x, −r ≤ θ ≤ 0. Definition 2 (Equilibria). The set of equilibria for Eq. (1) is defined by E {φ ∈ C+ : φ = xˆ for some x ∈ R+n satisfying f (x) ˆ = 0}. Definition 3 (Omega and Alpha Limit Sets). The omega limit set is defined by ω(φ) {x t (φ) : t ≥ s},
(3)
(4)
s≥0
whereas the alpha limit set is defined by α(φ) {x t (φ) : t ≤ s}.
(5)
s≤0
Definition 4 (Periodic Orbit). The orbit O+ (φ) is said to be a T-periodic orbit if x T +t (φ) = x t (φ) for all t and minimal T > 0. Next, we describe the cyclic feedback networks by mainly following the definitions and assumptions of Mallet-Paret and Sell (1996a,b), Kobayashi et al. (2003) and Wang et al. (2004a,b).
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
343
Assumption 2.1 (Monotone System). ∂ f i /∂ x j > 0, 0 (or 2, then this path is said to be a feedback loop. In addition, this feedback loop is said to be positive (or negative) if l−1 m=1 s pm+1 pm = 1 (or −1). In Fig. 2(a), there are many feedback loops, e.g., 1 → 2 → 1, 2 → 3 → 2 and 1 → 2 → 3 → 2 → 1, etc., but the largest loop, i.e., 1 → 2 → 3 → 4 → 5 → 1, is unique. Notice that Definition 8 requires l > 2. Actually, when l = 2, it is a self-feedback loop, and not included in the feedback loops in Definition 8 of the paper. Each node may have a linear or nonlinear, positive or negative self-feedback loop with no time delay in the network. All self-feedback loops are omitted in Figs. 1–3. Kobayashi et al. proved that a general biological network or Eq. (1) with only positive feedback loops has no dynamical attractors except stable equilibria (Kobayashi et al., 2002, 2003) when the system is bounded. In other words, except stable equilibria, there are neither stable periodic oscillations nor other nonequilibrium attractors for the network with only positive feedback loops. On the other hand, it is also indicated that a scalar negative delayed feedback system may have periodic solutions (Diekmann et al., 1995). According to Assumptions 2.1–2.3, Eq. (6) has only one largest feedback loop connecting all nodes, which indicates the nature of the whole system: positive or negative. Based on the monotone dynamical system theory and discrete Lyapunov functional, Mallet-Paret and Sell (1996a,b) obtained the Morse decomposition and derived the following theorem for the monotone cyclic feedback system (6) in the form of the Poincaré–Bendixson type theorem, when Assumptions 2.1–2.3 are satisfied. Let the natural phase space for Eq. (6) be C(K), where K = 0 τ21 ∪ 0 τ32 ∪ · · · ∪ 0 τn,n−1 ∪ 0 τ1,n ∪ N.
(11)
Theorem 2.1 (Poincaré–Bendixson Type Theorem). Assume that Assumptions 2.1–2.3 hold. Consider the system (6) satisfying Eq. (7) and being differentiable. Let x(t) be a solution of Eq. (6) on some time interval [t0 , ∞). Let ω(x) ⊆ C(K) denote the omegalimit set of this solution in the phase space C(K). Then either
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
347
Fig. 2. A simple case of the transformation procedure. All these three network structures have the same dynamical properties. (a) A feedback network with five nodes which satisfies Assumptions 2.1–2.4. (b) Types or signs of all interaction connected with the second node are changed. (c) On the basis of (b), types or signs of all interaction connected with the third node are changed. All self-feedback loops are omitted from the figures.
348
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Fig. 3. A scheme of a modified minimal model for circadian oscillations in PER and per mRNA with a negative feedback loop. Time delays are added into the original model to consider slow diffusion, transportation or signal transduction process of molecules between nucleus and cytosol.
1. ω(x) is a single nonconstant periodic orbit; or else 2. for each solution u(t) of Eq. (6) in ω(x), i.e., for solutions with u t ∈ ω(x) for all t ∈ R, we have α(u) ∪ ω(u) ⊆ E,
(12)
where α(u) and ω(u) denote the alpha- and omega-limit sets, respectively, of this solution, and where E ⊆ C(K) denotes the set of equilibria of Eq. (6). This theorem does not provide sufficient conditions for periodic orbits, but indicates that omega-limit sets of monotone cyclic feedback systems are composed of only periodic orbits and equilibria, which is a desirable property for modelling periodic oscillations of biological systems. In particular, the omega limit set is not empty due to the bounded assumption of Assumption 2.1. However, according to Assumption 2.3, the original networks are also restricted to special feedback structures, which may limit the potential applications. Motivated by Eq. (6) and Theorem 2.1, in this paper, we give sufficient conditions for a stable nonconstant periodic orbit and further extend our results to more general cyclic feedback networks, which have less restriction and are easy for applications. Eq. (6) or (1) with only positive feedback loops falls into the class of cooperative dynamical systems, which have been intensively investigated for both ODE and differential delay equations (Smith, 1991, 1995). Cooperative systems exhibit very regular behavior, e.g., typical solutions tend to equilibria in the case of autonomous systems, or to periodic solutions in the case of nonautonomous periodic systems. In other words, there are no stable periodic solutions for autonomous monotone dynamical systems with only positive
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
349
feedback loops (Smith, 1991, 1995; Gouze, 1998; Kobayashi et al., 2003). However, when there is a negative feedback loop in the monotone dynamical network, dynamics of Eq. (1) or (6) may become very complicated and possess a variety of nonequilibrium attractors, even in the autonomous cases. In particular, it can be shown that there is at most one equilibrium in X if all loops are nonpositive and the determinant of the Jacobian matrix of f is nonzero (Gouze, 1998). Therefore, to ensure the existence of the periodic solutions, we make the following important assumption. Assumption 2.4. The feedback loop connecting all nodes of the model for Eq. (6) is negative. Define the largest-loop to be the feedback loop connecting all nodes. Note that the largest-loop is unique for Eq. (6). When the largest-loop is negative, we say that Eq. (6) is a negative cyclic feedback network. In the paper, we prove that there exist attracting nonconstant periodic orbits for negative cyclic feedback networks, when the feedback for total one-direction interactions is sufficiently strong. In other words, we show that the omega-limit set is nonempty and is composed of only periodic orbits, which are quite different from those of cooperative systems. Fig. 2(a) is an example of an interaction graph with a negative largest-loop. The largest feedback loop, i.e., 1 → 2 → 3 → 4 → 5 → 1, is a negative feedback loop because of s21 s32 s43 s54 s15 = −1. According to Assumptions 2.2 and 2.3 and Fig. 1, clearly Theorem 2.1 for the negative cyclic feedback systems requires that 1. the largest-loop is negative and all self-feedback loops can have arbitrary signs; 2. all loops excluding the self-feedback loops and the largest-loop, are positive; 3. except the nth node, the interactions are opposite for the neighboring nodes, e.g., if e21 is +, then e32 must be −, as shown in Fig. 2(a). In this paper, besides sufficient conditions of the periodic orbits, we further extend Theorem 2.1 to general cyclic feedback systems, by eliminating 3 of the above conditions or Assumption 2.3. 3. Main results 3.1. Convergence to periodic solutions Cyclic nearest neighbor systems or cyclic feedback systems in the form of differential delay equations, i.e., Eq. (6), in which the coupling between neighbors possesses a monotonicity property, were extensively studied in Mallet-Paret and Sell (1996a,b). One important property of such systems is that the Poincaré–Bendixson type theorem holds for them. Therefore, an omega-limit set for such a system is either nonconstant periodic orbits or an equilibria if the system is bounded. When such systems are written in standard feedback forms, they are either positive (except self-feedback loops, all feedback loops including the largest-loop are positive) or negative (except self-feedback loops and negative largest-loop, all feedback loops are positive) feedback systems. In the case of positive feedback systems, such systems fall into the realm of the so-called cooperative systems.
350
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Orbits of cooperative systems or systems with all positive feedback loops have a strong tendency to converge to equilibria (Smith, 1995; Kobayashi et al., 2003). In contrast, for the negative feedback systems, attracting periodic orbits may exist, and their dynamics are quite different from cooperative systems. In this section, we derive our main theorems, and show that a cyclic feedback network can be used for modelling or designing a biological oscillator when the largest-loop is negative and the feedback for the total one-direction interactions is strong enough. Moreover, we further extend the original cyclic feedback networks to general cyclic feedback networks, which have less restriction. All detailed proofs for theoretical results are given in Appendix A. Since there is no limitation for the dimensionality, a bio-oscillator can be modelled or designed even by a large-scale system. Time delays do not change the location of an equilibrium in Eq. (6) but its stability. If all eigenvalues of the characteristic equation have negative real parts, then the equilibrium is stable and there is no oscillation locally. On the other hand, with change of a parameter, e.g., time delay τ , if any complex eigenvalue crosses the imaginary axis, then a stable equilibrium loses its stability with appearance of oscillation. Let x¯ ∈ X be an equilibrium of (6). Define
(13) for i, j ∈ N. Clearly, A(0) is the Jacobian matrix of f with respect where f i j = ∂∂xfij x=x¯ to x. Then the characteristic equation of Eq. (6) evaluated at the equilibrium x¯ is det(λI − A(λ)) = 0,
(14)
where I is the n × n identity matrix, and has the following form bn λn + bn−1 λn−1 + · · · + b0 + (−1)n+1 Be−λτ = 0, (15) n−1 where n is the number of nodes, B = f 1,n i=1 f i+1,i . b j for j = 0, . . . , n − 1 are functions of f k,k+1 f k+1,k and fi,i for 1 ≤ k ≤ n − 1 and i ∈ N except f 1,n , and bn = (−1)n . Total time delay is τ = i∈N τi+1,i . B represents the total feedback strength. Notice that b j for j = 0, . . . , n − 1 are functions of fk,k+1 f k+1,k for k ∈ N, which means that all effects of interactions between nodes k and k + 1 on b j disappear but on B exist if fk,k+1 is zero. Such a property is used in the proof of Theorem 3.2 in the appendix. From Eq. (15), clearly all time delays can be equivalently reduced to a single delay τ when stability of an equilibrium is considered, which is also indicated in Eq. (10) by the transformation of Eq. (8). Such a property drastically simplifies the analysis of Eq. (6) although there are multiple time delays. When the total time delay τ = 0, assume that the equilibrium x¯ is stable for Eq. (15). By changing the total time delay τ , we study the asymptotic stability of x¯ from the roots of the characteristic Eq. (15).
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
351
Since the solution exists and is bounded in X from the assumption, the omega limit set of Theorem 2.1 is nonempty, which means that all asymptotical solutions are periodic orbits provided that there is no stable equilibrium in X. Therefore, the basic idea in this paper is to destabilize all equilibria in X, which can be actually implemented mainly by linear analysis for most cases. We first state a theoretical result for the local analysis with respect to initial values, which guarantees a negative cyclic feedback system to converge to a local periodic orbit, by using the characteristic Eq. (15). Let mod (x, y) denote remainder after x divided by y. Theorem 3.1 (Local Convergence to Nontrivial Periodic Orbits). Assume that Assumptions 2.1–2.4 hold for Eq. (6), and the following equation has one nonzero real root v¯ at x¯ ∈ X, 2 2 i i−1 i i (−1) 2 bi v + (−1) 2 bi v − B 2 = 0, (16) i∈Ie
i∈Io
and further, 1 k−1 i−1 sgn (−1) 2 bk v¯ k (−1) 2 i bi v¯ i 2 2 B v¯ k∈Io i∈Io k i−2 − (−1) 2 bk v¯ k (−1) 2 i bi v¯ i = 0 k∈I i∈I ≥2 e
(17)
e,i
where Ie = {i ∈ N : mod (i, 2) = 0} ∪ {0}, Io = {i ∈ N : mod (i, 2) = 1} and x¯ is stable at τ = 0. Then there is a τ¯ , i (−1)n (−1) 2 bi v i 1 i∈Ie (18) τ¯ = arccos , v¯ B where the range of arccos is [0, π], and the system converges to a stable periodic orbit when τ is near τ¯ and τ > τ¯ . √ Theorem 3.1 implies that if there exists a τ¯ , at which, λ = jv¯ with j = −1 is a root of Eq. (15) and the derivative of the real part of the eigenvalue at λ = jv¯ is not zero, then the cyclic feedback system will converge to a stable periodic orbit when τ is near τ¯ and τ > τ¯ , for any initial conditions near x. ¯ This theorem shows that a biological network with a negative feedback loop can have periodic orbits, which can be proven by the Hopf bifurcation theorem for functional differential equations (Kolmanovskii and Myshkis, 1999). Moreover, we show that the Hopf bifurcation is supercritical due to the stable equilibrium x¯ before bifurcation, as indicated in Appendix A. Notice that there is no restriction on the self-feedback loop for cyclic feedback systems, i.e., for any i ∈ N, the interaction from the i th chemical component directly to itself may be either nonlinear or linear, and either positive or negative. Although the expressions of Theorem 3.1 seem to be complex, the conditions are easy to satisfy.
352
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Based on the local convergence of Theorem 3.1, we describe the global convergence conditions of existence for nontrivial periodic orbits. Theorem 3.2 (Global Convergence to Nontrivial Periodic Orbits). If Eq. (6) satisfies Assumptions 2.1–2.4 and the feedback for total one-direction interactions is sufficiently ∂f ∂ fi strong, i.e., ∂∂xf1n i ∂i+1 x i for those i with ∂ x i+1 = 0 at any equilibrium is sufficiently large, then there exists a τ¯ , i (−1)n (−1) 2 bi v i 1 i∈Ie (19) τ¯ = arccos , v¯ B where the range of arccos is [0, π], and the system converges to a stable periodic orbit for almost all initial conditions φ ∈ Xˆ , when τ > τ¯ . ∂ f 1 ∂ f i+1 ∂ fi i ∂ x i for all those i with ∂ x i+1 = 0 represents the total strength of one-direction ∂ xn interactions, i.e., the product of those interactions ∂ f i+1 /∂ x i with si,i+1 = 0 but si+1,i = 0. Theorem 3.2 is proven by showing that all equilibria are unstable due to the positive real part of an eigenvalue and the real part of that eigenvalue never returns to a negative region for any τ > τ¯ . Therefore, there is no stable equilibrium but stably periodic orbits in X due to the nonempty omega limit set of Eq. (6). Actually, we can prove that when p = 2 or 3, only b02 − B 2 < 0 is needed for existence of stable periodic orbits in Theorem 3.2, instead of sufficiently strong total one-direction interactions. Theorem 3.3 (Global Convergence to Nontrivial Periodic Orbits). Assume that Eq. (6) satisfies Assumptions 2.1–2.4. If det(A(0)) < 0 at all equilibria in X, then for almost all initial conditions φ ∈ Xˆ , x t (φ) converges to a stable periodic solution, where A(λ) is defined by (13). The conditions of Theorem 3.3 are strong, and ensure that the system converges to a stable periodic orbit, regardless of any nonnegative time delays. Theorem 3.3 is proven by showing that all equilibria in X are unstable for any nonnegative time delays. It is generally not easy to guarantee stable behaviors such as equilibria or periodic orbits even for a small network with a few components, due to the nonlinearity of the system. Theorem 3.2 means that if the feedback for total one-direction interactions is sufficiently strong, an attracting periodic orbit exists. In other words, the omega-limit set is nonempty and composed of only periodic orbits. Moreover, according to the proof of Theorem 3.2 in Appendix A, the total time delay τ affects the dynamics of Eq. (6), in contrast to the cooperative systems. On the other hand, Theorem 3.3 indicates that when the determinant of the Jacobian matrix is negative for all equilibria in X, the system converges to periodic orbits with any initial function φ ∈ Xˆ and any nonnegative time delays. As indicated in Kobayashi et al. (2003), if all feedback loops of a system are positive, then its orbits have a strong tendency to converge to equilibria. However, for negative cyclic feedback systems, when conditions of Theorem 3.2 or 3.3 hold, attracting periodic orbits exist and only periodic orbits constitute the omega-limit sets, which are quite different
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
353
from those of positive feedback networks. Therefore, if conditions of Theorem 3.2 or 3.3 are satisfied, Eq. (6) will inevitably converge to stable periodic orbits. In other words, negative cyclic feedback systems have ideal properties for oscillating networks and can be used to model and design bio-oscillators. Although Theorem 3.1 is a local convergence theorem, it becomes a global convergence result for any initial condition when there is at most one equilibrium in X, as stated in Corollary 3.1. Corollary 3.1 (Global Convergence to Nontrivial Periodic Orbits). Assume that all conditions of Theorem 3.1 hold. If det(A(0)) = 0 for all x in a convex set X, then the system converges to a stable periodic orbit for almost all initial condition φ ∈ Xˆ , when τ is near τ¯ and τ > τ¯ . By showing that there is at most one equilibrium that is unstable, we can prove the corollary, which is given in the appendix. 3.2. Generalized cyclic feedback systems Generally sign (feedback) conditions on the nonlinearities f i in Assumption 2.3 require that interactions are opposite for the neighboring nodes except the first and last nodes, which are difficult to satisfy for many biological systems. Next, we eliminate such restrictions by coordinate transformations. In other words, we extend the cyclic feedback systems to general ones with any type of interactions for the neighboring nodes. For Eq. (6), which satisfies Assumptions 2.1–2.3, choose any node i ∈ N and change types or signs of all interactions connected with the node-i , we denote the system obtained under such a transformation as y˙ (t) = g(yt )
(20)
where the transformation P for x = Py is defined by a matrix as follows: for the node i∈N 0 σ1 .. (21) P= , . 0
σn
with σi = −1 and σ j = 1 for all j but j = i . By substituting x = Py into Eq. (20), we have another FDE as follows: x(t) ˙ = f (x t ) ≡ Pg(Px t )
(22)
P−1
where P = is used. We can easily prove that Eq. (22) satisfies all Assumptions 2.1–2.3. Moreover, since P is a reversible and one-to-one map, Eq. (22) is qualitatively equivalent to Eq. (20). Therefore, we have the following theorem. Theorem 3.4 (General Cyclic Feedback Systems). For Eq. (6), assume that Assumptions 2.1–2.3, hold. The transformation of Eq. (21), which changes signs of all interactions connected with any node i (1 ≤ i ≤ n), does not change its dynamical properties. Theorem 3.4 means that Eq. (20), in which the signs of all interactions connected with any node i ∈ N are changed, is qualitatively equivalent to Eq. (22), which satisfies
354
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Assumptions 2.1–2.3. Moreover, it is also easy to show that such transformation does not change the type of any feedback loop, which implies that a negative cycle feedback network is still a negative cycle feedback network under this transformation. A simple case of the transformation procedure is shown in Fig. 2(a)–(c). The difference between the Fig. 2(a) and (b) is types or signs of all interactions connected with the second node. All cases in Fig. 2 are dynamically equivalent. An important property of this transformation is that we can get any combination of interactions under such consecutive transformations, which actually do not change type of each loop. By performing such transformations for each node, we can obtain all combinations of cyclic feedback networks with any type of interactions, which are all qualitatively equivalent. Therefore, Assumption 2.3 can be eliminated for the general cyclic feedback networks, by the following corollary. Corollary 3.2. If Assumptions 2.1 and 2.2 are satisfied, then Theorems 2.1 and 3.1–3.3 and Corollary 3.1 still hold for Eq. (6). In contrast to the special structures of the original cyclic feedback networks for Theorem 2.1, Corollary 3.2 means that general negative cyclic feedback systems with any type of interactions on each node can be used to design oscillators provided that the largestloop is negative. Fig. 2 shows different feedback network structures, which have the same dynamical properties. The results of general cyclic feedback networks actually are further applied to a hybrid monotone dynamical system with both positive and negative feedback loops (Wang et al., 2004b). 4. Numerical implementation 4.1. A minimal model for circadian rhythm with time delays In this section, we demonstrate our theoretical results by modifying a minimal model (Goldbeter, 1995) for circadian rhythms in the Drosophila, which uses a period protein (PER) and the period per gene. In particular, we add time delays into the original model (Goldbeter, 1995) to consider slow diffusion, transportation or signal transduction process of molecules between nucleus and cytosol. per mRNA, whose cytosolic concentration is denoted by M, is assumed to be synthesized in the nucleus and transferred to the cytosol with time delay τm for the translation of PER, where it is degraded. The rate of synthesis of PER is proportional to M. Three states of the protein are considered in this model, i.e., unphosphorylated (P0 ), mono-phosphorylated (P1 ) and bisphosphorylated (P2 ). In this model, we consider that PN which is a complex of P2 , behaves directly as a repressor. per mRNA affected from PN with time delay τ p , is synthesized in the nucleus and transferred to the cytosol, where it accumulates at a maximum rate vs ; and it is degraded by an enzyme with maximum rate vm and Michaelis constant K m . The rate of synthesis of the PER protein from M with time delay τm , is characterized by a first-order rate constant ks . Parameters Vi and K i (i = 1, . . . , 4) denote the maximum rate and Michaelis constant of the kinase(s) and phosphatase(s) involved in the reversible phosphorylation of P0 into P1 and P1 into P2 , respectively. The repressor PN is synthesized and degraded by the first-order rate constant k2 and k1 in cytoplasm, and
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
355
is transported into the nucleus to inhibit the mRNA in the Hill type with time delay τ p . The scheme of the model for circadian oscillations in PER and per mRNA is shown in Fig. 3. The time evolution of the five-variable model is governed by the following functional differential equations, in which all parameters and concentrations are defined with respect to the total cell volume, P1 (t) P0 (t) + ks M(t − τm ) − V1 P˙0 (t) = V2 K 2 + P1 (t) K 1 + P0 (t) P0 (t) P2 (t) P˙1 (t) = V1 + V4 − D1 (P1 (t)) K 1 + P0 (t) K 4 + P2 (t) P1 (t) (23) + k2 PN (t) − D2 (P2 (t)) P˙2 (t) = V3 K 3 + P1 (t) P˙N (t) = k1 P2 (t) − k2 PN (t) K In M(t) ˙ − vm M(t) = vs n K I + PNn (t − τ p ) K m + M(t) where D2 (P2 (t)) = V4
P2 (t) P2 (t) + k1 P2 (t) + vd K 4 + P2 (t) K d + P2 (t)
D1 (P1 (t)) = V2
P1 (t) P1 (t) + V3 . K 2 + P1 (t) K 3 + P1 (t)
and
The total quantity of PER protein, Pt is given by Pt = P0 + P1 + P2 + PN ,
(24)
where n denotes the degree of cooperativity, and K I is the threshold constant for repression. By using the transformation to a single time delay, i.e., Eqs. (8) and (23) can be equivalently changed to a system with only one time delay with identical forms as the original ones except that M(t − τm ) → M (t)
(25)
and PN (t − τ p ) is replaced by PN (t − τ ) in the last equation of Eq. (23), where total time delay is τ = τ p + τm . In the following discussion, we will use the new variable M in (25). All parameters are from Goldbeter (1995) with slight modifications, and are set as vs = 0.9125 µM h−1 , vm = 1.625 µM h−1 , K m = 0.5 µM, ks = 2.25 h−1 , vd = 2.5 µM h−1 , k1 = 6.5 h−1 , k2 = 1 h−1 , K I = 1 µM, K d = 0.1 µM, n = 4, K 1 = 0.65, K 2 = 2, K 3 = 1, K 4 = 2 µM, V2 = 3.95 µM h−1 , V3 = 12.5 µM h−1 , V4 = 6.25 µM h−1 , V1 = 10. τ (µM h−1 ) is given in each case next. Clearly Eq. (23) is a negative cyclic feedback network and satisfies Assumptions 2.1 and 2.2. Fig. 4 shows a case for a sustained oscillation generated by Eq. (23) with τ = 8, for which the period of the oscillation is around 24 h. Fig. 5 shows the trajectory toward a limit cycle as a projection onto the plane formed by the concentrations of per mRNA M and total PER protein Pt . Because perturbations do not change the period or amplitude in
356
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Fig. 4. A sustained oscillation generated from Eq. (23) based on the negative feedback regulation of per mRNA by the PER protein in Drosophila at τ = 8. The period of the oscillation is around 24 h.
Fig. 5. Evolution toward a limit cycle corresponding to a sustained oscillation shown in Fig. 4 at τ = 8. The arrow indicates the direction of movement along the orbit.
the long run, a limit cycle oscillation represents a particularly stable mode of the periodic behavior. Such stability holds with the robust nature of circadian clocks which have to maintain their amplitude and period in the changing environment. According to Eqs. (16) and (18) of Theorem 3.1, we can obtain v¯ = 0.3949 and τ¯ = 4.78, which are bifurcation values and also consistent with our numerical results.
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
357
Fig. 6. A bifurcation diagram with total time delay τ as a parameter. The solid line and dashed lines represent the equilibrium is stable or unstable respectively. The dash-dotted lines indicate the maximum and minimum values of M for the sustained oscillations. Limit cycles exist when τ > τ¯ , for which the equilibrium in this region is unstable.
Since Theorem 3.1 gives local bifurcation conditions for a periodic orbit, a sustained oscillation occurs after a bifurcation point by some control parameters. This situation is shown in Fig. 6, where the control parameter is the total time delay τ . At low values of τ , the system reaches a stable steady state. With the increase of τ , a bifurcation occurs at a critical τ¯ = 4.78. After τ > τ¯ , the steady state becomes unstable and a sustained oscillation appears. The amplitudes of the sustained oscillations are also shown in Fig. 6. When τ = τ¯ , we get a pair of imaginary roots λ = ±0.3949j for the characteristic equation. Moreover, from Fig. 6, the amplitudes increase with the time delay τ , which means that the time delay affects the amplitudes of the oscillations. When τ < τ¯ , the system reaches a stable steady state. For example, Fig. 7 shows that a trajectory asymptotically converges to an equilibrium in the (Pt , M ) plane at τ = 4.70 < τ¯ . Fig. 8 is the analysis of the effect for total time delay τ on the period T . In addition to amplitudes as shown in Fig. 6, the period of oscillation, namely T , increases with the total time delay τ in almost a linear way. Therefore, the total time delay τ can be viewed as a key parameter to control both amplitude and period of an oscillation in a biological system. 4.2. Repressilator with time delays We next examine a synthetic gene regulatory network, i.e., the Repressilator by genes cI, TetR and LacI, whose periodic oscillation was experimentally investigated in a prokaryote, i.e., E. coli in vivo (Elowitz and Leibler, 2000). For eukaryotes, due to slow processes of transportation, diffusion or signal transduction of molecules between nucleus and cytoplasm, the time delays may play an important role for the cellular dynamics. In this paper, we assume that the Repressilator is constructed in an eukaryote, such as yeast. The model can be expressed as (Elowitz and Leibler, 2000; Chen and Aihara, 2002a; Chen et al., 2004):
358
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Fig. 7. An asymptotical trajectory to a stable equilibrium in the (Pt , M ) plane at τ = 4.7.
Fig. 8. Total time delay τ and period of oscillations T .
α − α1 m i (t) 1 + p j (t − τ p j )n p˙ i (t) = βm i (t − τm i ) − βpi (t)
m˙ i (t) =
(26) (27)
where i and j have the following three pairs of values: (i = 1, j = 2), (i = 2, j = 3) and (i = 3, j = 1). n is the Hill coefficient. m i ∈ R+ are mRNAs, and pi ∈ R+ are proteins. τ p j and τm i are time delays. The regulation scheme of each gene is illustrated in Fig. 9.
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
359
Fig. 9. Repressilator with three mRNAs m i and three proteins pi , (i = 1, 2, 3), signs + and − on the edges mean s = 1 and −1, respectively. Self-feedback loops are omitted from the figure.
Let total time delay be τ=
3 i=1
τm i +
3
τ pi .
(28)
i=1
By the transformation Eq. (8), that is m 1 (t − τ + τ p2 ) → m 1 (t),
m 2 (t − τm 2 ) → m 2 (t), m 3 (t − τm 3 − τ p3 − τm 2 ) → m 3 (t), p1 (t − τ + τ p2 + τm 1 ) → p1 (t), p2 (t) → p2 (t),
(29)
p3 (t − τ p3 − τm 2 ) → p3 (t),
Eqs. (26) and (27) are equivalently converted into the following form α − α1 m i (t), m˙ i (t) = 1 + p j (t − τ p j )n p˙ i (t) = βm i (t) − βpi (t),
(30) (31)
for (i = 1, j = 2), (i = 2, j = 3) and (i = 3, j = 1), where τ p2 = τ and τ p1 = τ p3 = 0. Parameters are set as n = 2, β = 2, α1 = 1 and α = 2.1. Clearly, Eqs. (26) and (27) satisfy Assumptions 2.1 and 2.2. From Eqs. (16) and (18) of Theorem 3.1, we obtain the bifurcation values v¯ = 0.1969 and τ¯ = 11.5, which are also confirmed by numerical simulation. The bifurcation diagram and analysis of the effect for total time delay τ on the period T are shown in Figs. 10 and 11 respectively, which obviously have the similar features as the circadian model. It is easy to check that b0 = 8 and |B| = 8.5975, which means that b02 − B 2 < 0. Therefore, the conditions of Theorem 3.2 are also satisfied, which
360
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Fig. 10. Effects of τ on amplitude of m 1 . The solid line and dashed lines represent stable and unstable equilibria respectively. The dash-dotted lines indicate the maximum and minimum values of m 1 for the sustained oscillations.
Fig. 11. Total time delay τ and period of oscillations T .
implies that there is no stable equilibria but stable periodic solutions for Eqs. (26) and (27) when the feedback loop is sufficiently strong. 5. Conclusion This paper developed a new methodology to analyze and design biological oscillating networks with time delays, by using cyclic feedback systems. One desirable property for a monotone cyclic feedback system is that there exists a nonempty omega-limit set composed
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
361
of only periodic orbits and equilibria. Such a property drastically simplifies theoretical analysis and design of oscillators. We first provided sufficient conditions for existence of limit cycles from both local and global analysis, and then further extended the original cyclic feedback systems into general cyclic feedback networks, which allow any types of interactions on each node. Such an extension relaxes the restrictions of cyclic feedback systems, and enables our model to apply to a wide variety of systems for modelling and designing bio-oscillators. Finally, two examples are used for numerical simulation, one from a circadian network (Goldbeter, 1995) formed by a period protein (PER) and per mRNA, and the other from a synthetic gene regulatory network, i.e., Repressilator (Elowitz and Leibler, 2000). Since time delays may play an important role particularly in eukaryotes, we added time delays into the models to study their effects on the cellular dynamics. We showed that a periodic oscillation is generally produced by the negative feedback in this model. Since the time delays may significantly influence the dynamics of the overall system, we also analyzed the effects of the total delay on the stability of the equilibria and periodic oscillations. Simulation results indicated that the total time delay can slowly change the amplitudes and vary periods in almost a linear manner. Our simulation results imply that total time delay can be used as an important control parameter both quantitatively and qualitatively to change dynamics of the systems, in particular for the amplitudes and periods of periodic oscillations. Although we have mainly examined effects of time delays on the cellular dynamics, there are also other important facts, which may play crucial roles in biological processes and should be further investigated in future works from both theoretical and experimental viewpoints, such as stochastic noise (Elowitz et al., 2002; Swain et al., 2002; Zhou et al., submitted for publication) and discrete nature of biological models. Acknowledgements The work is partly supported by Osaka Sangyo University, and National Natural Science Foundation of China (10371037). Appendix A. Proofs of theorems and corollary In this appendix, we give detailed proofs of all theorems and corollary in this paper. A.1. Proof of Theorem 3.1 With the Assumptions 2.1–2.4, Eq. (1) can be written as Eq. (6), which is a monotone negative cyclic feedback system with delays. The roots of the transcendental Eq. (15) for λ determine the stability of equilibria. If real parts of all roots are negative for an equilibrium, then it is stable. On the other hand, if there exists a root with a positive real part, the equilibrium becomes unstable. √ In other words, a Hopf bifurcation occurs when λ = jv is a root of Eq. (15), where j = −1 and v
362
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
is a nonzero real number, that is i (−1) 2 bi v i = (−1)n cos(vτ )B, i∈Ie
(−1)
i−1 2
(A.1)
bi v i = (−1)n+1 sin(vτ )B,
i∈Io
where Ie = {i ∈ N : mod (i, 2) = 0} ∪ {0} and Io = {i ∈ N : mod (i, 2) = 1}. For any τ , by eliminating terms of sin and cos, we have
2 i 2
(−1) bi v
i
+
i∈Ie
2 (−1)
i−1 2
bi v
− B 2 = 0.
i
(A.2)
i∈Io
Hence, if Eq. (A.2) has a nonzero real root v, ¯ then Eq. (15) has an imaginary root λ = jv. ¯ Notice that B < 0 due to the negative cyclic feedback loop. Therefore, we can obtain the critical values for v¯ and τ¯ , which has the following form τ¯ =
1 arccos v¯
(−1)n
i
(−1) 2 bi v¯ i
i∈Ie
B
,
(A.3)
where the range of arccos is [0, π]. On the other hand, from Eq. (15), we have
dλ dτ
−1
= (−1)n+1
h (λ) λτ τ e − , Bλ λ
(A.4)
where h (λ) = (−1)n nλn−1 + (n − 1)bn−1 λn−2 + · · · + b1 .
(A.5)
Note that i−1 i−2 1 i i (−1) 2 i bi v¯ + j (−1) 2 i bi v¯ h (λ)|λ=jv¯ = v¯ i∈I i∈I ,i≥2
o
1 (h 1 + j h 2 ). v¯
e
(A.6)
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
363
λ) dλ −1 Since sgn{ d(Re dτ } = sgn Re{( dτ ) }, by substituting Eq. (A.1) into the real part of Eq. (A.4) at λ = jv, ¯ we obtain ∂(Re λ) sgn ∂τ λ=jv¯ n+1 h (λ) λτ e = sgn Re (−1) Bλ λ=jv¯ ! h sin( vτ ¯ ) + h ¯ ) 2 cos(vτ n+1 1 = sgn (−1) (A.7) B v¯ 2 k−1 i−1 1 k 2 bk v = sgn (−1) ¯ (−1) 2 i bi v¯ i B 2v¯ 2 k∈I i∈Io o k i−2 k i . − (−1) 2 bk v¯ (−1) 2 i bi v¯ k∈Ie
i∈Ie ,i≥2
Thus, when the right part of Eq. (A.7) is nonzero, we have
∂(Re λ) ∂τ λ=jv¯
> 0, which implies
that real parts of the complex eigenvalues move to the right half plane for increasing τ with λ on the imaginary axis. Then the stable equilibrium x¯ loses its stability, according to Hopf bifurcation theorem for delay equations. Therefore, Eq. (6) has a stable periodic solution when τ is near τ¯ and τ > τ¯ . In other words, the Hopf bifurcation is supercritical. In fact, by assumption, x¯ is a stable equilibrium at τ = 0, which means when 0 < τ < τ¯ , λ) < 0 means there all real parts of eigenvalues are negative. On the other hand, ∂(Re ∂τ λ=jv¯
are at least a pair of complex eigenvalues with positive real parts when 0 < τ < τ¯ , which contradicts our assumption. This completes the proof. A.2. Proof of Theorem 3.2 Assume that feedback for total one-direction interactions is sufficiently strong, i.e., ∂ fi+1 ∂ fi i ∂ x i for those i with ∂ x i+1 = 0 at any equilibrium is sufficiently large. Then, according to the definition of B from Eq. (15), |B| is also sufficiently large, which implies that ∂ f1 ∂ xn
b02 − B 2 < 0.
(A.8)
In other words, at any equilibrium, the even-order Eq. (A.2) at least has a nonzero real root v, ¯ which means that characteristic Eq. (15) has an imaginary root λ = jv. ¯ Note that B < 0 due to the assumption of the negative cyclic feedback network. On the other hand, if B > 0, Eq. (6) is actually a network with only positive feedback loops, which has no stable periodic orbit in X according to Kobayashi et al. (2003) and Smith (1995). Moreover, using Eqs. (A.1) and (A.2), we have the following relationship for sufficiently large |B|: 1
v¯ ∼ |B| n ,
(A.9)
364
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
and for k = 0, 1, 2, . . .. τ¯k =
1 2kπ + arccos v¯
(−1)n
i
(−1) 2 bi v¯ i
i∈Ie
B
,
where the range of arccos is [0, π]. From Eq. (A.4), we have −1 dλ h (λ) λτ τ e − , = (−1)n+1 dτ Bλ λ and d(Re λ) sgn dτ
! = sgn Re
dλ dτ
−1
(A.10)
(A.11)
.
Therefore, we obtain n v¯ 2n−2 d(Re λ) = sgn sgn (1 + O( )) , dτ λ=jv¯ B2
(A.12)
(A.13)
where the expression O( ) = O( n1 ((b12 −2b0b2 )v¯ 2−2n +c2n−4 v¯ 4−2n +· · ·+c4 v¯ −4 +c2 v¯ −2 )) represents a higher order of , which is small when B is sufficiently large, where ck (k = 2, 4, . . . , 2n − 4) are functions of bi , (i = 0, 1, . . . , n − 1). When the feedback for total one-direction interactions is sufficiently strong, |B| is d(Re λ) > 0 for any number of nodes. Therefore, sufficiently large. Thus we have dτ λ=jv¯
the real part of λ moves to the right half plane for τ > τ¯ when λ is on the imaginary axis, i.e., the stable equilibrium loses the stability. However, by the assumption, the solution of Eq. (1) or (6) exists and is bounded in X, which implies that there exists a nonempty omega limit set composed of only periodic orbits and equilibria according to the Poincaré–Bendixson theorem for cyclically coupled differential delay equations (Mallet-Paret and Sell, 1996a,b) or Theorem 2.1. Since such arguments hold for all equilibria, all equilibria lose their stability and the omega limit sets are composed of only stable periodic orbits, which proves this theorem. A.3. Proof of Theorem 3.3 Consider the characteristic Eq. (15) evaluated at an equilibrium x¯ det(λI − A(λ)) = 0,
(A.14)
where I is the n × n identity matrix, and has the following form h(λ) ≡ bn λn + bn−1 λn−1 + · · · + b1 λ + b0 + (−1)n+1 Be−λτ = 0,
(A.15)
where n is the number of nodes and bn = (−1)n . Assume A(0) < 0. We have h(0) = b0 − B = A(0) < 0 and h(+∞) = +∞ for an even number n; h(0) = b0 + B = −A(0) > 0 and h(+∞) = −∞ for an odd number n. By continuality of h(λ), h(λ) = 0 has at least
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
365
one positive real root for both of the above situations. Therefore, x¯ is unstable. If the above conditions hold at all equilibria, then all equilibria are unstable. By assumption, the solution of Eq. (6) exists and is bounded in X. Therefore, according to the Poincaré–Bendixson theorem for systems of cyclically coupled differential delay equations (Mallet-Paret and Sell, 1996a,b) or Theorem 2.1, there is a nonempty omega limit set composed of only stable periodic orbits, which proves this theorem. A.4. Proof of Theorem 3.4 For Eq. (6), which satisfies Assumptions 2.1–2.3, choose any node i ∈ N and change types or signs of all interactions connected with the node-i . We denote the system obtained under such a transformation as y˙ (t) = g(yt )
(A.16)
where the transformation P is defined by a matrix as follows: for the node i ∈ N 0 σ1 .. P= , . 0
(A.17)
σn
with σi = −1 and σ j = 1 for all j but j = i . By substituting x = Py into Eq. (A.16), we have another FDE as follows: x(t) ˙ = f (x t ) = Pg(Px t )
(A.18)
where P = P−1 is used. According to Eq. (A.17), the mapping from g to f is one to one and reversible. Therefore, f is topologically equivalent to g and Eq. (A.18) has the same dynamical properties as Eq. (A.16). Moreover, because of ∂ f i /∂ x i−1 = σi σi−1 ∂gi /∂yi−1 with σi = −1 and σi−1 = 1, ∂ f i /∂ x i−1 and ∂gi /∂yi−1 have opposite signs. In the same way, it is easy to check that ∂ f i+1 /∂ x i , ∂ f i−1 /∂ x i and ∂ f i /∂ x i+1 also change their signs. However, all other feedbacks, do not change their signs. Clearly Eq. (A.18) satisfies all our Assumptions 2.1 and 2.2. This completes the proof. A.5. Proof of Corollary 3.1 To prove the corollary, we only need to show that f is injective on the convex set X or there is at most one equilibrium, when A(0) is not singular for all x ∈ X according to Gouze (1998). Assume that x 1 and x 2 are two different points in the convex set X such that f (x 1 ) = f (x 2 ). Then, " 1 D f ((1 − α)x 1 + αx 2 )(x 2 − x 1 )dα 0 = f (x 1 ) − f (x 2 ) = 0
"
1
=
A(0)dα(x 2 − x 1 )
0
where D f = A(0) is the Jacobian matrix of f with respect to x.
(A.19)
366
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
Clearly, if A(0) is not singular for all x ∈ X, then the determinant det(A(0)) is nonzero and has a constant sign due to continuity of the system. Since (1 − α)x 1 + αx 2 is also #1 in the convex set X for 0 ≤ α ≤ 1, det( 0 A(0)dα) is nonzero or the integrated matrix #1 0 A(0)dα is nonsingular, which implies x 1 = x 2 from Eq. (A.19). This contradicts the assumption, and proves the corollary. References Barbai, N., Leibler, S., 2000. Biological rhythms: circadian clocks limited by noise. Nature 403, 267. Becskei, A., Serrano, L., 2000. Engineering stability in gene networks by autoregulation. Nature 405, 590–593. Belair, J., Campbell, S.A., Driessche, P., 1996. Frustration, stability, and delay-induced oscillations in a neural network model. SIAM J. Appl. Math. 56, 245–255. Campbell, S.A., Ruan, S., Wei, J., 1999. Qualitative analysis of a neural network model with multiple time delays. Int. J. Bifurcation Chaos 9, 1585–1595. Chen, L., Aihara, K., 2002a. Stability of genetic regulatory networks with time delay. IEEE Trans. Circuits Syst. I. 49, 602–608. Chen, L., Aihara, K., 2002b. A model of periodic oscillation for genetic regulatory systems. IEEE Trans. Circuits Syst. I. 49, 1429–1436. Chen, L., Wang, R., Kobayashi, T., Aihara, K., 2004. Dynamics of gene regulatory networks with cell division cycles. Phys. Rev. E 70, 011909, 1–13. Crosthwaite, S.K., Dunlap, J.C., Loros, J.J., 1997. Neurospora wc-1 and wc-2: transcription, photoresponses, and the origin of circadian rhythmicity. Science 276, 763–769. Diekmann, O., van Gils, S.A., Verduyn Lunel, S.M., Walther, H.O., 1995. Delay Equations Functional-, Complex-, and Nonlinear Analysis. Springer-Verlag. Dunlap, J.C., 1999. Molecular bases for circadian clocks. Cell 96, 271–290. Elowitz, M.B., Leibler, S., 2000. A synthetic oscillatory network of transcriptional regulators. Nature 403, 335–338. Elowitz, M.B., Levine, A.J., Siggia, E.D., Swain, P.S., 2002. Stochastic gene expression in a single cell. Science 297, 1183–1186. Gardner, T.S., Cantor, C.R., Collins, J.J., 2000. Construction of a genetic toggle switch in Escherichia Coli. Nature 403, 339–342. Glass, L., Kauffman, S.A., 1973. The logical analysis of continuous, non-linear biochemical control networks. J. Theor. Biol. 39, 103–129. Goldbeter, A., 1995. A model for circadian oscillations in the Drosophila period protein (PER). Proc. R. Soc. Lond. B 261, 319–324. Golden, S.S., Ishiura, M., Johnson, C.H., Kondo, T., 1997. Cyanobacterial circadian rhythms. Ann. Rev. Plant Physiol. Plant Mol. Biol. 48, 327–354. Gonze, D., Leloup, J.-C., Goldbeter, A., 2000. Theoretical models for circadian rhythms in Neurospora and Drosophila. C. R. Acad. Sci. Paris, III 323, 57–67. Goodwin, B.C., 1965. Oscillatory behavior in enzymatic control process. Adv. Enzyme Regua. 3, 425–438. Gouze, J.-L., 1998. Positive and negative circuits in dynamical systems. J. Biol. Syst. 6, 11–15. Hall, J.C., 1998. Genetics of biological rhythms in Drosophila. Adv. Genet. 38, 135–184. Hasty, J., Isaacs, F., Dolnik, M., McMillen, D., Colins, J.J., 2001. Designer gene networks: towards fundamental cellular control. Chaos 11, 207–220. Hunding, A., 1974. Limit cycles in enzyme systems with nonlinear negative feedback. Biophys. Struct. Mech. 1, 47–54. Kobayashi, T., Chen, L., Aihara, K., 2002. Design of genetic switches with only positive feedback loops. In: Proceedings of 2002 IEEE Computer Society Bioinformatics Conference. pp. 151–162. Kobayashi, T., Chen, L., Aihara, K., 2003. Modelling genetic switches with positive feedback loops. J. Theor. Biol. 221, 379–399. Kolmanovskii, V., Myshkis, A., 1999. Introduction to the Theory and Applications of Functional Differential Equations. Kluwer Academic Publishers.
R. Wang et al. / Bulletin of Mathematical Biology 67 (2005) 339–367
367
Leloup, J.-C., Goldbeter, A., 1998. A model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins. J. Biol. Rhythms 13, 70–87. Leloup, J.-C., Goldbeter, A., 2000. Modelling the molecular regulatory mechanism of circadian rhythms in Drosophila. BioEssays 22, 84–93. Mallet-Paret, J., Sell, G.R., 1996a. Systems of differential delay equations: floquet multipliers and discrete Lyapunov Functions. J. Differ. Eqns. 125, 385–440. Mallet-Paret, J., Sell, G.R., 1996b. The Poincaré–Bendixson theorem for monotone cyclic feedback systems with delay. J. Differ. Eqns. 125, 441–489. Mikkar, A.J., Carre, I.A., Strayer, C.A., Chua, N.H., Kay, S.A., 1995. Circadian clock mutants in Arabidopsis identified by iuciferase imaging. Science 267, 1161–1163. Smith, H., 1991. Periodic tridiagonal competitive and cooperative systems of differential equations. SIAM J. Math. Anal. 22, 1102–1109. Smith, H., 1995. Monotone Dynamical Systems, vol. 41. American Mathematical Society, Providence, RI. Swain, P.S., Elowitz, M.B., Siggia, E.D., 2002. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc. Natl. Acad. Sci. 99, 12795–12800. Tei, H., Okamura, H., Shigeyoshi, Y., Fukuhara, C., Ozawa, R., Hirose, M., Sakaki, Y., 1997. Circadian oscillation of a mammalian homologue of the Drosophila period gene. Nature 389, 512–516. Walther, H.O., 1995. The 2-dimensional attractor of x(t) = −µx(t) + f (x(t − 1)). Memoirs of the Amer. Math. Soc., 544. Amer. Math. Soc, Providence, RI. Wang, R., Jing, Z., Chen, L., 2004a. Periodic oscillators in genetic networks with negative feedback loops. WSEAS Trans. Math. 3, 150–156. Wang, R., Zhou, T., Jing, Z., Chen, L., 2004b. Modelling periodic oscillation of biological systems with multiple time scale networks. Syst. Biol. 1, 71–84. Wu, J., Zou, X., 1995. Patterns of sustained oscillations in neural networks with delayed interactions. Appl. Math. Comput. 73, 55–75. Zhou, T., Chen, L., Aihara, K., 2004. Intercellular communications induced by random fluctuations (submitted for publication).
