- Email: [email protected]
Discrete and ultradiscrete models for biological rhythms comprising a simple negative feedback loop
Journal of Theoretical Biology 378 (2015) 89–95
Contents lists available at ScienceDirect
Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi
Discrete and ultradiscrete models for biological rhythms comprising a simple negative feedback loop Shingo Gibo a,nn,1, Hiroshi Ito b,n a b
Graduate School of Design, Kyushu University, 4-9-1, Shiobaru Minami-ku, Fukuoka 815-8540, Japan Faculty of Design, Kyushu University, 4-9-1, Shiobaru Minami-ku, Fukuoka 815-8540, Japan
H I G H L I G H T S
We propose discrete model with two variables comprising negative feedback. Discrete model generates self-sustained oscillations for some conditions. We derive ultradiscrete model from discrete model. Ultradiscrete model generates self-sustained oscillations for some conditions. Ultradiscrete model includes a Boolean system as a special case.
art ic l e i nf o
a b s t r a c t
Article history: Received 28 December 2014 Received in revised form 14 April 2015 Accepted 20 April 2015 Available online 30 April 2015
Many biological rhythms are generated by negative feedback regulation. Griffith (1968) proved that a negative feedback model with two variables expressed by ordinary differential equations do not generate self-sustained oscillations. Kurosawa et al. (2002) expanded Griffith's result to the general type of negative feedback model with two variables. In this paper, we propose discrete and ultradiscrete feedback models with two variables that exhibit self-sustained oscillations. To obtain the model, we applied tropical discretization and ultradiscretization to a continuous model with two variables and then investigated its bifurcation structures and the conditions of parameters for oscillations. We found that when the degradation rate of the variables is lower than their synthesis rate, the proposed models generate oscillations by Neimark–Sacker bifurcation. We further demonstrate that the ultradiscrete model can be reduced to a Boolean system under some conditions. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Self-sustained oscillations Biological rhythms Neimark–Sacker bifurcation Boolean system Tropical discretization
1. Introduction Numerous examples of self-sustained oscillations are known at the molecular to ecological levels, including cell division cycles, segmentation clocks (Hirata et al., 2002), oscillations of p53 expression (Lahav et al., 2004), and periodical outbreaks of insects (Liebhold et al., 2000). Many of these biological rhythms are generated by negative feedback loops. At the molecular level, autonegative feedback loops of clock genes, such as per, frq, cca1 and kai, produce circadian rhythms (Aronson et al., 1994; Ishiura et al., 1998; McWatters and Devlin, 2011; Partch et al., 2014). Negative feedback loops in a genetic circuit are very common such
n
Corresponding author. Tel./fax: þ 81 925534535. Corresponding author. Tel.: +81 484679546; fax: +81 484621709 E-mail addresses: [email protected] (S. Gibo), [email protected] (H. Ito). 1 Present address: Theoretical Biology Laboratory, RIKEN, 2-1 Hirosawa, Wako 351-0198, Japan nn
http://dx.doi.org/10.1016/j.jtbi.2015.04.024 0022-5193/& 2015 Elsevier Ltd. All rights reserved.
that they are well known as a network motif (Alon, 2007). More than 40% transcription factors in Escherichia coli repress their own gene expression (Shen-Orr et al., 2002; Rosenfeld et al., 2002). At the macroscopic level, a predator–prey system in an ecosystem can be regarded as a negative feedback loop. The oscillations in population size have been intensively examined from not only theoretical but also experimental viewpoints (Fussmann et al., 2000). The range of number of reactions that constitute negative feedback loops is broad. Some repressors can directly repress their own expression. The HES1 protein, which is responsible for the segmentation clock, can directly repress its own expression by binding to the promoter of the hes1 gene (Hirata et al., 2002). In contrast, the circadian clock genes form a complex negative feedback loop with several steps, including multiple phosphorylation and nuclear transportation (Partch et al., 2014). In this paper, we focus on the negative feedback loop comprising two elements in biological networks. Griffith (1968) proved that no oscillations are generated by the negative feedback model
90
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
a _ ¼ described as M bM; E_ ¼ cM dE through Bendixson 1 þ KEm criterion. Kurosawa et al. (2002) expanded Griffith's model to the generalized case. They proved that no oscillations are produced by the general class of negative feedback loop described as follows: dx1 ¼ f ðx2 Þ hðx1 Þ; dt
ð1aÞ
dx2 ¼ gðx1 Þ kðx2 Þ; dt
ð1bÞ
where f ðx2 Þ is a monotonically decreasing function, and gðx1 Þ, hðx1 Þ, and kðx2 Þ are monotonically increasing functions. The proof was given by the construction of the Lyapunov function for Eqs. (1). This model corresponds to the case of a common feedback loop of a single gene through transcription and translation. Suppose x1 to be the amount of messenger RNA and x2 to be that df of protein. is negative, given that proteins decrease the rate of dx2 dg is positive, given that gene expression via negative feedback. dx1 dh 4 0 and proteins are produced via the translation of mRNA. dx1 dk 4 0 should hold because of the degradation of the mRNA and dx2 proteins. These conditions should also be satisfied in the simple case of a predator–prey relationship. When x1 and x2 represent the df populations of prey and predator, respectively, o 0 and dx2 dg dh dk 4 0 hold owing to predation, and 4 0 and 4 0 hold dx1 dx1 dx2 because of their death. In the negative feedback with n variables as follows: dx1 ¼ f 1 ðxn Þ h1 ðx1 Þ; dt dx2 ¼ f 2 ðx1 Þ h2 ðx2 Þ; dt
ð2aÞ
ð2bÞ
⋮ dxn ¼ f n ðxn 1 Þ hn ðxn Þ; dt the necessary condition for instability is β 1 β 2 ⋯β n π n 4 sec ; α 1 α 2 ⋯α n n
ð2cÞ
ð3Þ
dhi ðxi Þ df ðx Þ ði ¼ 1; 2; …; nÞ, β i ¼ i i 1 ði ¼ 2; 3; …; nÞ and dxi dxi 1 df ðx Þ β1 ¼ 1 n . This condition is known as the secant condition dxn (Tyson and Othmer, 1978; Thron, 1991; Sontag, 2006). In particular,
where αi ¼
if n ¼2, any values of αi and βi cannot satisfy the condition (3) because the right-hand side of Eq. (3) is infinity. This result is consistent with Kurosawa et al. (2002). The proof that negative feedback with two variables cannot generate oscillation is limited to ordinary differential equations (ODEs). Some researchers have recently focused on cases where ODEs are not applicable. For example, gene expression noise is unavoidable, owing to the low copy numbers of most genes. The effect of molecular noise on the oscillation was assessed via stochastic simulation (Barkai and Leibler, 2000; Gonze et al., 2002; Nishino et al., 2013). Although impossibility of oscillations was proved for ODE model of negative feedback with two variables, whether the other model for negative feedback with
two variables can show self-sustained oscillations has not yet been investigated. Besides ODEs, discrete and ultradiscrete systems are used to model dynamical systems. A discrete system treats the independent variables such as time and space as discrete numbers and is represented as a difference equation (Elaydi, 2005). An ultradiscrete system treats both independent and dependent variables as discrete numbers (Tokihiro et al., 1996). The qualitative behaviors may not generally be conserved by discretizing or ultradiscretizing ODEs. For example, although the logistic equation described as an ODE is attracted to an equilibrium point, the logistic map obtained by the forward discretization of the logistic equation shows a chaotic behavior (May, 1976). In the present study, we focus on discrete and ultradiscrete negative feedback models with two variables. These models can oscillate under some parameter conditions. We also investigated the types of bifurcation and parameter conditions for the oscillations.
2. Discretization In this section, we discretize a negative feedback model with two variables and study its bifurcation structure. 2.1. Model We consider the negative feedback model with two variables, which is reduced from the three-variable model proposed by Kurosawa et al. (2002), as follows: dM ¼ dt
k n aM; P 1þ h
ð4aÞ
dP ¼ sM vP: dt
ð4bÞ
If Eqs. (4) represent a genetic feedback loop, then M and P correspond to amounts of mRNA and protein, respectively. We assume that the parameters k, a, s, v, n, and h are positive. Because Eqs. (4) are included in the class of (1), Eqs. (4) have no limit cycle solutions. To nondimensionalize Eqs. (4), we define the dimenk ks 1 sionless variables M ¼ x, P ¼ 2 y, t ¼ τ and dimensionless v v v 2 a hv parameters b ¼ and c ¼ . Then, Eqs. (4) are transformed into v ks dx 1 ¼ ð5aÞ y n bx; dτ 1þ c dy ¼ x y: dτ
ð5bÞ
By applying tropical discretization (Murata, 2013) to Eqs. (5), we obtain 1 y n bxT 1þ T xT þ 1 xT c ¼ ; δ 1 þ bδ yT þ 1 yT
δ
¼
xT yT ; 1þδ
ð6aÞ ð6bÞ
or, xT þ xT þ 1 ¼
δ
1þ
y n T
c
1 þ bδ
;
ð7aÞ
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
91
Fig. 1. Numerical simulations on the discrete system (Eqs. (7)); (a) n¼2, b¼ 1, c ¼ 0.01 and δ ¼ 3, (b) n¼2, b¼ 1, c ¼0.01 and δ ¼ 1.
yT þ 1 ¼
yT þ δxT ; 1þδ
ð7bÞ
where δ is the time interval. We see that this discretization is relevant because if we substitute τ ¼ T δ and take the continuous limit δ-0, we obtain Eqs. (5). The numerical simulation confirms that the discrete model (7) appears to oscillate or to converge to an equilibrium point depending on parameters n, b, c, and δ (Fig. 1).
From Eq. (11) and b; δ 40, the parameters of n and b, and xn must satisfy n
n 1 4 nbx :
ð12Þ
Note that when δ is sufficiently large, the discrete model (7) has a limit cycle as long as n, b, and xn satisfy Eq. (12). Given that n, b, and xn are positive, parameter n must satisfy n 4 1:
2.2. Stability
We next consider the oscillatory parameter region of b and c.
The parameter region for oscillation can be calculated by stability analysis. Eqs. (7) have a unique fixed point, ðxn ; yn Þ, which satisfies n
bx ¼
ð13Þ
1 n n ; y 1þ c
ð8aÞ
n
Combining Eqs. (8a) and (8b), we obtain bx þ
n
ðbx Þn þ 1 ¼ 1; this is ðbcÞn
rewritten as bc ¼
" #1=n n ðbx Þn þ 1 : n 1 bx
ð14Þ n
xn ¼ yn :
ð8bÞ
Since 0 o
1þ
1yn n o 1, the range of value of xn is c
1 0 o xn o : b
ð9Þ
The Jacobian at ðxn ; yn Þ is 2 3 n 1 nbδ1ð1þbδbx Þ 1 þ bδ 4 5: A¼ δ 1 1þδ
ð10Þ
1þδ
þ bÞδ and The trace and the determinant of A are τ ¼ ð12þþbð1δÞð1 þ δÞ þ nbδ ð1 bx Þ Δ ¼ 1 ð1 , þ bδÞð1 þ δÞ n
2
pffiffiffiffiffiffiffiffiffiffiffiffi
respectively,
and
the
eigenvalues
are
λ ¼ τ 7 τ2 4Δ. 2
If τ2 4Δ o 0, the eigenvalues are complex. τ 2 4Δ is negative for a broad range of values of b and δ, as can be confirmed numerically (Appendix A). In particular, if b is nearly unity, then τ2 4Δ is always negative. We now consider the case of τ2 4Δ o 0. pffiffiffiffi When j λ j ¼ Δ 4 1, i.e., ½n 1 nbx bδ 4 1 þ b; n
From the oscillatory condition of Eq. (12), we obtain bx o n n 1. Combining them, the parameters b and c under oscillatory conditions must satisfy
ð11Þ
the fixed point loses its stability, and the discrete model (7) generates self-sustained oscillations via Neimark–Sacker bifurcation. We can derive that τ2 4Δ o 0 is always satisfied when Eq. (11) holds (Appendix B); i.e., self-sustained oscillation occurs via Neimark–Sacker bifurcation whenever the fixed point is destabilized.
bc o
ðn 1Þðn þ 1Þ=n : n
ð15Þ
Thus, when parameters n, b, c, and δ satisfy Eqs. (11), (13), and (15), the fixed point ðxn ; yn Þ is unstable, and the discrete system (7) behaves as an oscillator. We next consider the range of the value of bc. Eq. (15) can be transformed into bc on1=n ð1 1nÞ1 þ ð1=nÞ . Moreover, n1=n ð1 1nÞ1 þ ð1=nÞ-1 as n-1. Hence, the value of bc is bounded because of the continuity of
ðn 1Þðn þ 1Þ=n . n
ðn 1Þðn þ 1Þ=n is below n 2 hv c ¼ ks , parameters
av h o 1:2: ks
By numerical computation, the maximum of
1.2. Recalling that b and c are defined by b ¼ av and a, v, k, s, and h satisfy ð16Þ
Recall that a and v denote the degradation rates of M and P, respectively, k and s are the synthesis rates of M and P, and h is the threshold for negative regulation. Thus, Eq. (15) suggests that selfsustained oscillation occurs when the degradation rate is lower than the synthesis rate and the threshold of feedback regulation is small. Note that the continuous limit of the Jacobian, limδ-0 A δ I, where I is the identity matrix, equals the Jacobian of the differential equations (5). This result means that if δ is sufficiently small, the period of the discrete system (7) equals that of the continuous system (5).
92
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
ultradiscretization method to Eqs. (22). Let us introduce
2.3. Generalization We will find a general condition for oscillation in discrete negative feedback model with two variables, regardless its function form. We consider the general form of the negative feedback with two variables, that is, Eqs. (1) substituted x1 ¼ x and x2 ¼ y as follows: dx ¼ f ðyÞ hðxÞ; dt
ð17aÞ
dy ¼ gðxÞ kðyÞ; dt
ð17bÞ dk 4 0, dy
where and respectively. By applying tropical discretization to Eqs. (17), we obtain the discrete model xT þ δf ðyT Þ ; xT þ δhðxT Þ
ð18aÞ
y þ δgðxT Þ yT þ 1 ¼ yT T : yT þ δkðxT Þ
ð18bÞ
The Jacobian at the fixed point ðxn ; yn Þ is 2 3 δ xn α 1 δ xn β 1 n 6 1 xn þ δhðxn Þ x þ δhðxn Þ 7 6 7 ; A¼6 n δy β 2 δyn α2 7 4 5 1 n n n n y þ δkðy Þ y þ δkðy Þ
ð19Þ
b¼e
;
yT ¼ eY T =ε ;
c ¼ eC=ε ;
δ ¼ eD=ε :
Then, we obtain
X T þ 1 ¼ X T þ ε log eX T =ε þ eD=ε
δ2 xn yn β1 β2 ; ðxn þ δhðxn ÞÞðyn þ δkðyn ÞÞ
ð24aÞ
X T þ 1 ¼ X T þ max½X T ; max½0; nðY T CÞ þ D max½X T ; X T þ B þ D;
ð25aÞ
Y T þ 1 ¼ Y T þ max½Y T ; X T þ D max½Y T ; Y T þD;
ð25bÞ
where we use the formula for ultradiscretization as follows: lim ε log ½eA1 =ε þeA2 =ε þ ⋯ þ eAn =ε ¼ max½A1 ; A2 ; …; An :
ε- þ 0
ð20Þ
Simulating the ultradiscretized model (25), we can observe oscillations when n Z1, D 4 max½ B; 0, and B þ C o 0 (Fig. 2). These conditions correspond to the ultradiscretization of the oscillatory conditions (11) and (15) for the discrete model, as shown below. From Eqs. (11) and (23), we obtain B=ε
We next focus on the ultradiscretization of our developed discrete model. Eqs. (7) can be rewritten as
yT þ 1 ¼ yT
δ
y n T
yT þ δxT : yT þ δyT
;
ð27Þ
D 4 max½ B; 0: Similarly, from Eqs. (15) and (23), we obtain " # ðn 1Þðn þ 1Þ=n B þ C o ε log : n
ð28Þ
ð29Þ
Taking the limit ε- þ 0, we obtain the ultradiscretized condition B þ C o 0:
ð30Þ
4. Reduction to a Boolean system
3. Ultradiscretization
xT þ 1 ¼ xT
þ BÞ=ε
þ 1 εlog ½n 1:
Assuming that the time interval is large, we can make the discrete model (7) and ultradiscrete model (25) simpler. Taking the limit δ-1 for the discrete model (7), we obtain 1 xT þ 1 ¼ y n ; b 1þ T c yT þ 1 ¼ xT ;
c xT þ δbxT
n
ð21Þ
If the time interval δ is nonzero value, the negative feedback model with two variables (18) might oscillate when the synthesis rates are larger than the degradation rates. By taking continuous limit δ-0, the right-hand side of Eq. (21) diverges to infinity and any values of α1, α2, β1 and β2 cannot satisfy Eq. (21). This means the continuous negative feedback model (17) cannot the generate oscillations. This result corresponds to the secant condition when n ¼2 (Tyson and Othmer, 1978; Thron, 1991; Sontag, 2006).
1þ
ð26Þ
Taking the limit ε- þ 0, we have the ultradiscretized condition
dh dk and α2 ¼ are the degradation rates of x and y, dx dy df dg respectively, β1 ¼ and β 2 ¼ are the synthesis rates of x and y. dy dx If Δ 4 1, the fixed point is unstable. Therefore, the sufficient condition for instability is
β1 β2 xn þ δhðxn Þ yn þ δkðyn Þ 4 þ 1: α1 α2 δα1 xn δα2 yn
ð24bÞ
Taking the limit ε- þ 0, we have the ultradiscrete equations
4 ε log ½e
where α1 ¼
xT þ
ð23Þ
1 ε log ½eX T =ε þ eX T þ B þ D=ε ; 1 þ eðnðY T CÞÞ=ε
D 4 ε log ½e B=ε þ 1 εlog ½n 1 neðX
and the determinant of A is δ xn α 1 δyn α2 1 zw Δ ¼ 1 n x þ δhðxn Þ yn þ δkðyn Þ þ
B=ε
Y T þ 1 ¼ Y T þ ε log ½eY T =ε þ eðX T þ DÞ=ε ε log ½eY T =ε þ eðY T þ DÞ=ε :
df o 0, dh 4 0, dg 40 dy dx dx
xT þ 1 ¼ xT
xT ¼ eX T =ε ;
ð22aÞ
ð22bÞ
Tokihiro et al. (1996) proposed a method for deriving ultradiscrete equations from discrete equations. We will apply the
ð31aÞ
ð31bÞ
where the right-hand side of Eq. (31a) represents a negative feedback effect. According to the simulation of Eqs. (31), this model can display oscillation. Thus, Eqs. (31) represent the simplest discrete model for self-sustained oscillations produced by a negative feedback loop. Similarly, taking the limit D-1 for the ultradiscrete model (25) X T þ 1 ¼ max½0; nðY T CÞ B;
ð32aÞ
Y T þ 1 ¼ XT ;
ð32bÞ
where the right-hand side of Eq. (32a) represents a negative feedback effect. The numerical simulation of Eqs. (32) confirms that this model, as well as the discrete system, conserves the oscillatory solution. The conditions of δ-1 and D-1 can be interpreted as a sufficiently long time interval between successive reactions in
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
93
Fig. 2. Numerical simulations on the ultradiscrete system (Eqs. (25)); (a) n¼ 1, B ¼0, C ¼ 1 and D ¼1, (b) n ¼1, B¼ 1, C ¼ 1 and D¼ 1.
biochemical systems or as cycles between long generation times in ecosystems. These simple models generally correspond to negative feedback loops where the reaction interval is relatively larger than reaction rate. Moreover, Eqs. (32) can be reduced to a Boolean system under the condition shown below. A Boolean system is one that consists of binary variables and logical operations, such as AND ( 4 ), OR ( 3 ), and NOT (:). Suppose n ¼1, C r 0, and C B ¼ 1. Eqs. (32) become X T ¼ :Y T ;
ð33aÞ
Y T ¼ XT ;
ð33bÞ
where 1 Y T is replaced with :Y T . This model is a Boolean system because the range of XT and YT is restricted in f0; 1g if we assign 0 or 1 to the variables as initial values. This means that the ultradiscrete model for negative feedback loops, Eqs. (25), includes a Boolean system, Eqs. (33), as a special case.
5. Discussion and conclusions We focused on the simple negative feedback model that consists of two elements. As Tyson and Othmer (1978), Thron (1991) and Kurosawa et al. (2002) proved, the ODE model with two variables never oscillates, regardless of its function form. We found that the discretized and ultradiscretized models can oscillate depending on the values of the parameters. In discrete model, we derived the general condition for oscillation, which covers any forms of function. Discretizing and ultradiscretizing can contribute to destabilizing a stable negative feedback system. We have also found the developed ultradiscrete model can be reduced to a Boolean system under some conditions. Some textbooks illustrate the mechanism of circadian rhythms produced by a simple negative feedback loop as follows: (i) expression of mRNA induces the production of proteins, (ii) proteins inhibit the expression via feedback, (iii) lowered expression level decreases protein amounts, (iv) the expression starts again with the disappearance of the inhibitory proteins (Alberts et al., 2008). This model is equivalent to the Boolean system Eqs. (33). Note that a discretized or ultradiscretized simple feedback model can show oscillations but that a continuous limit of the model cannot always conserve its rhythmicity. This observation may explain why most biological rhythms consist of many molecular parts. The parameter n represents the nonlinearity of negative feedback control. n corresponds to an index of the cooperativity of protein–DNA interaction in the context of a gene regulatory network. Some studies have shown that higher cooperativity is more likely to contribute in the destabilization of the negative
feedback system described in an ODE (Leloup and Goldbeter, 1998; Kurosawa et al., 2002). For example, the negative feedback model with three variables proposed in Kurosawa et al. (2002) can be destabilized when n Z 9. In contrast, our developed model can be destabilized when n 4 1 for the discrete model and n Z 1 for the ultradiscrete model, suggesting that the nonlinearity of feedback is not necessary for discretized oscillators. We used the tropical discretization for discretizing Eqs. (5). This discretization was proposed by Murata (2013) as a systematic manner of discretizing and ultradiscretizing first-order differential equations. This discretization preserves the positivity of variables if the initial values and parameters are positive. Moreover, dependent variables x and y are proportional to the amounts of some proteins or the populations of prey and predator. Thus, x and y must be positive, and tropical discretization is appropriate in this case. In contrast, applying the other discretization to Eqs. (5), the system may not preserve positivity. For example, applying the forward difference to Eqs. (5), the discrete model has a negative value when δ is large. We consider the biological interpretation of the tropical discretization. The right-hand sides of Eqs. (6) are the right-hand sides of Eqs. (5) divided by 1 þ bδ or 1 þ δ. Thus, the larger the time-step size δ, the smaller are the right-hand sides of Eqs. (6). This result means that when the time interval between a chemical reaction and the next chemical reaction is large, the effect of synthesis and the degradation of a product are small. This assumption may be valid for real biological systems to preserve the positivity of components. As mentioned in the Introduction, the stochastic model can be an alternative to reproduce the behavior of dynamical systems in nature. We transformed the deterministic model (4) into a stochastic model (Appendix C). The obtained model produced pulses intermittently for a broad range of the values of parameters. This occurs because the state drifting around the fixed point happens to be ejected due to the molecular noise and randomly makes a large-amplitude rotation. Unlike the discrete and ultradiscrete model, we did not succeed in observing any oscillations with a constant period and amplitude. The discretized or ultradiscretized model corresponds to some actual biological systems. Behavior as a discretized system is often essential to an ecological system. Many species of plants and animals perform reproductive activity every year. The population dynamics can be described as a discrete system with δ ¼1 [year]. It is the common case that the production rate of species Y is enhanced by species X and conversely the species Y decrease the production rate of species X. Based on our research, we can conclude that this simple ecological negative feedback regulation can cause the oscillation as long as the condition (21) is satisfied. The number of population should be represented as a discrete value rather than continuous one when the distribution of individuals is sparse and
94
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
Appendix A. Condition for τ2 4Δ o 0
12
We numerically computed the sign of τ2 4Δ, which is the discriminant of the characteristic polynomial of the Jacobian shown in Eq. (10) (Fig. 3). We found that τ2 4Δ is negative for most values of b; δ. In particular, τ2 4Δ is always negative when b is near b¼1, resulting in self-sustained oscillation or damped oscillation. Recall that b is defined as b ¼ av and that a and v are degradation rates of M and P, respectively. Hence, b ¼1 means that the degradation rate of M nearly equals that of P.
10 8 6 4
Appendix B. The condition for τ2 4Δ o 0 includes the condition for Δ 41
2 0
0
2
4
6
8
10
12
Numerical calculation suggests that Δ 4 1 is satisfied when τ2 4Δ o 0 holds and when n ¼ 3; c ¼ 0:8 (Fig. 3). Moreover, the area of τ2 4Δ o0 includes the area of Δ 4 1 for arbitrary n and δ. This is analytically confirmed as follows:
Number of
Fig. 3. Sign of τ2 4Δ. The black area represents the sets of b and c that satisfy τ2 4Δ 40 when n¼ 3 and c¼ 0.8. The broken line denotes Δ ¼ 1 at which Neimark– Sacker bifurcation occurs.
τ2 4Δ ¼ ½ð2 þ ð1 þ bÞδÞ2 4f1 þ nbδ2 ð1 bxn Þgð1 þ bδÞð1 þ δÞ =ð1 þbδÞ2 ð1 þ δÞ2
200
¼ ½ 4nb ð1 bx Þδ 4nbð1 þbÞð1 bx Þδ
100
Assuming Δ 41, we obtain δ 4
0
n
2
n
2
4nbð1 bx Þ þ ð1 bÞ δ =ð1 þ bδÞ ð1 þ δÞ2 : n
0
20
40
60
Time Fig. 4. Stochastic simulations on the system Eqs. (25); h ¼ a ¼ v ¼ 1, k ¼ v ¼ 5, Ω¼ 100 and τ¼ 1.
ðB:1Þ
2
2
2
ðB:2Þ
1þb . By substituting the n ½n 1 nbx b inequality into the second term of the numerator in Eq. (B.2), we have 1þb τ2 4Δ o 4nb2 ð1 bxn Þδ2 4nbð1 þ bÞð1 bxn Þ n ½n 1 nbx b i 2 n 4nbð1 bx Þ þ ð1 bÞ2 δ =ð1 þ bδÞ2 ð1 þ δÞ2 ðB:3Þ "
ð1 bÞ2 n n 1 nbx n
nð3b þ 1Þðb þ 3Þð1 bx Þ 2 δ =ð1 þ bδÞ2 ð1 þ δÞ2 : þ n n 1 nbx
¼ 4nbð1 bx Þðbδ þ 1Þ þ n
the number of species is quite low. Ultradiscrete method should be applicable in this case. At the molecular level, biochemical reaction can be also described using the discrete model due to time delay. For example, intron produces the time delay in a transcriptional process. It is theoretically well known that the time delay can contribute to destabilize the system. Some experimental research confirms this hypothesis. Hes7 gene that is responsible for segmentation clock forms autonegative feedback loop, and shows the self-sustained oscillations of gene expression. The intronic region of Hes7 is required for the genetic oscillations (Takashima et al., 2011; Harima et al., 2013). In addition, the length of intron can lengthen the period of transcriptional oscillation in engineered genetic oscillator (Swinburne et al., 2008). Meanwhile, post-translational modification and nuclear transportation of proteins can generate time delay in circadian rhythms. The discrete model can be used to explicitly incorporate the time delay between synthesis and being mature molecule that can function. The ultradiscrete method should be used when the number of molecules is low as well as the ecological level. The effect of a low number of individuals or molecules needs to be addressed through an ultradiscrete model in future studies.
n
2
ðB:4Þ
n
Because 1 bx 4 0 and n 1 nbx 4 0 can be derived from Eqs. (9) and (12), respectively, we conclude τ2 4Δ o 0 if Δ 4 1. Appendix C. Stochastic simulation for genetic feedback loop with two variables We checked the rhythmicity of the stochastic model transformed from the deterministic feedback model with two variables. If Eqs. (4) represent a genetic feedback loop, the ODE model can be decomposed into some elementary biochemical reactions that consist of a gene G, mRNA M and protein P as follows: G þ nP-GP n ;
ðC:1aÞ
GP n -G þ nP;
ðC:1bÞ
G-G þM;
ðC:1cÞ
M-∅;
ðC:1dÞ
M-M þ P;
ðC:1eÞ
Acknowledgment
P-∅:
ðC:1f Þ
We would like to thank T. Kawabe (Kyushu University) for many helpful discussions. This research was supported by Inamori Foundation and Takeda Science Foundation (H.I.).
The two tional to time of (Nishino
scale parameters Ω and τ are introduced; Ω is proporthe number of molecules in the system, and τ scales the association and dissociation of proteins to the DNA et al., 2013). We assumed the transition rates of the
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
reactions (C.1a)–(C.1f) are 1τ ðPhΩÞn G, 1τ ð1 GÞ, kΩG, aM, sM and vP. We performed stochastic simulation for the reactions using the Gillespie algorithm and observed stochastic pulse for a broad range of parameters (Fig. 4).
References Alberts, B., Johnson, A., Lewis, J., Raff, M., Keith, R., Peter, W., 2008. Molecular Biology of the Cell, 5th ed. Garland Science, New York. Alon, U., 2007. Network motifs: theory and experimental approaches. Nat. Rev. Genet. 8 (6), 450–461. Aronson, B.D., Johnson, K.A., Loros, J.J., Dunlap, J.C., 1994. Negative feedback defining a circadian clock: autoregulation of the clock gene frequency. Science 263, 1578–1584. Barkai, N., Leibler, S., 2000. Biological rhythms: circadian clocks limited by noise. Nature 403, 267–268. Elaydi, S., 2005. An Introduction to Difference Equations, 3rd ed. Springer, New York. Fussmann, G.F., Ellner, S.P., Shertzer, K.W., Hairston, N.G., 2000. Crossing the Hopf bifurcation in a live predator–prey system. Science 290, 1358–1360. Gonze, D., Halloy, J., Goldbeter, A., 2002. Robustness of circadian rhythms with respect to molecular noise. Proc. Natl. Acad. Sci. USA 99, 673–678. Griffith, J.S., 1968. Mathematics of cellular control processes I. Negative feedback to one gene. J. Theor. Biol. 20, 202–208. Harima, Y., Takashima, Y., Ueda, Y., Ohtsuka, T., Kageyama, R., 2013. Accelerating the tempo of the segmentation clock by reducing the number of introns in the Hes7 gene. Cell Rep. 3, 1–7. Hirata, H., Yoshiura, S., Ohtsuka, T., Bessho, Y., Harada, T., Yoshikawa, K., Kageyama, R., 2002. Oscillatory expression of the bHLH factor Hes1 regulated by a negative feedback loop. Science 298, 840–843. Ishiura, M., Katsuna, S., Aoki, S., Iwasaki, H., Andersson, C.R., Tanabe, A., Golden, S.S., Johnson, C.H., Kondo, T., 1998. Expression of a gene cluster kaiABC as a circadian feedback process in cyanobacteria. Science 281, 1519–1523. Kurosawa, G., Mochizuki, A., Iwasa, Y., 2002. Comparative study of circadian clock models, in search of processes promoting oscillation. J. Theor. Biol. 216, 193–208.
95
Lahav, G., Rosenfeld, N., Sigal, A., Geva-Zatorsky, N., Levine, A.J., Elowitz, M.B., Alon, U., 2004. Dynamics of the p53-Mdm2 feedback loop in individual cells. Nat. Genet. 36, 147–150. Leloup, J.C., Goldbeter, A., 1998. Modeling circadian oscillations of PER and TIM proteins in Drosophila. In: Touitou, Y. (Ed.), Biological Clocks. Mechanisms and Applications. Elsevier Science, Amsterdam, pp. 81–88. Liebhold, A., Elkinton, J., Williams, D., Muzika, R.M., 2000. What causes outbreaks of the gypsy moth in North America? Popul. Ecol. 42, 257–266. May, R.M., 1976. Simple mathematical models with very complicated dynamics. Nature 261, 459–467. McWatters, H.G., Devlin, P.F., 2011. Timing in plants—a rhythmic arrangement. FEBS Lett. 585 (10), 1474–1484. Murata, M., 2013. Tropical discretization: ultradiscrete Fisher-KPP equation and ultradiscrete Allen–Cahn equation. J. Differ. Equ. Appl. 19 (6), 1008–1021. Nishino, R., Sakaue, T., Nakanishi, H., 2013. Transcription fluctuation effects on biochemical oscillations. PLoS ONE 8 (4), e60938. Partch, C.L., Green, C.B., Takahashi, J.S., 2014. Molecular architecture of the mammalian circadian clock. Trends Cell Biol. 24 (2), 90–99. Rosenfeld, N., Elowitz, M.B., Alon, U., 2002. Negative autoregulation speeds the response times of transcription networks. J. Mol. Biol. 323 (5), 785–793. Shen-Orr, S.S., Milo, R., Mangan, S., Alon, U., 2002. Network motifs in the transcriptional regulation network of Escherichia coli. Nat. Genet. 31, 64–68. Sontag, E.D., 2006. Passivity gains and the “secant condition” for stability. Syst. Control Lett. 55, 177–183. Swinburne, I.A., Miguez, D.G., Landgraf, D., Silver, P.A., 2008. Intron length increases oscillatory periods of gene expression in animal cells. Genes Dev. 22, 2342–2346. Takashima, Y., Ohtsuka, T., González, A., Miyachi, H., Kageyama, R., 2011. Intronic delay is essential for oscillatory expression in the segmentation clock. Proc. Natl. Acad. Sci. USA 108 (8), 3300–3305. Thron, C.D., 1991. The secant condition for instability in biochemical feedback control—I. The role of cooperativity and saturability. Bull. Math. Biol. 53 (3), 383–401. Tokihiro, T., Takahashi, D., Matsukidaira, J., Satsuma, J., 1996. From soliton equations to integrable cellular automata through a limiting procedure. Phys. Rev. Lett. 76 (18), 3247–3250. Tyson, J.J., Othmer, H.G., 1978. The dynamics of feedback control circuits in biochemical pathways. Prog. Theor. Biol. 5, 1–63.
Contents lists available at ScienceDirect
Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi
Discrete and ultradiscrete models for biological rhythms comprising a simple negative feedback loop Shingo Gibo a,nn,1, Hiroshi Ito b,n a b
Graduate School of Design, Kyushu University, 4-9-1, Shiobaru Minami-ku, Fukuoka 815-8540, Japan Faculty of Design, Kyushu University, 4-9-1, Shiobaru Minami-ku, Fukuoka 815-8540, Japan
H I G H L I G H T S
We propose discrete model with two variables comprising negative feedback. Discrete model generates self-sustained oscillations for some conditions. We derive ultradiscrete model from discrete model. Ultradiscrete model generates self-sustained oscillations for some conditions. Ultradiscrete model includes a Boolean system as a special case.
art ic l e i nf o
a b s t r a c t
Article history: Received 28 December 2014 Received in revised form 14 April 2015 Accepted 20 April 2015 Available online 30 April 2015
Many biological rhythms are generated by negative feedback regulation. Griffith (1968) proved that a negative feedback model with two variables expressed by ordinary differential equations do not generate self-sustained oscillations. Kurosawa et al. (2002) expanded Griffith's result to the general type of negative feedback model with two variables. In this paper, we propose discrete and ultradiscrete feedback models with two variables that exhibit self-sustained oscillations. To obtain the model, we applied tropical discretization and ultradiscretization to a continuous model with two variables and then investigated its bifurcation structures and the conditions of parameters for oscillations. We found that when the degradation rate of the variables is lower than their synthesis rate, the proposed models generate oscillations by Neimark–Sacker bifurcation. We further demonstrate that the ultradiscrete model can be reduced to a Boolean system under some conditions. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Self-sustained oscillations Biological rhythms Neimark–Sacker bifurcation Boolean system Tropical discretization
1. Introduction Numerous examples of self-sustained oscillations are known at the molecular to ecological levels, including cell division cycles, segmentation clocks (Hirata et al., 2002), oscillations of p53 expression (Lahav et al., 2004), and periodical outbreaks of insects (Liebhold et al., 2000). Many of these biological rhythms are generated by negative feedback loops. At the molecular level, autonegative feedback loops of clock genes, such as per, frq, cca1 and kai, produce circadian rhythms (Aronson et al., 1994; Ishiura et al., 1998; McWatters and Devlin, 2011; Partch et al., 2014). Negative feedback loops in a genetic circuit are very common such
n
Corresponding author. Tel./fax: þ 81 925534535. Corresponding author. Tel.: +81 484679546; fax: +81 484621709 E-mail addresses: [email protected] (S. Gibo), [email protected] (H. Ito). 1 Present address: Theoretical Biology Laboratory, RIKEN, 2-1 Hirosawa, Wako 351-0198, Japan nn
http://dx.doi.org/10.1016/j.jtbi.2015.04.024 0022-5193/& 2015 Elsevier Ltd. All rights reserved.
that they are well known as a network motif (Alon, 2007). More than 40% transcription factors in Escherichia coli repress their own gene expression (Shen-Orr et al., 2002; Rosenfeld et al., 2002). At the macroscopic level, a predator–prey system in an ecosystem can be regarded as a negative feedback loop. The oscillations in population size have been intensively examined from not only theoretical but also experimental viewpoints (Fussmann et al., 2000). The range of number of reactions that constitute negative feedback loops is broad. Some repressors can directly repress their own expression. The HES1 protein, which is responsible for the segmentation clock, can directly repress its own expression by binding to the promoter of the hes1 gene (Hirata et al., 2002). In contrast, the circadian clock genes form a complex negative feedback loop with several steps, including multiple phosphorylation and nuclear transportation (Partch et al., 2014). In this paper, we focus on the negative feedback loop comprising two elements in biological networks. Griffith (1968) proved that no oscillations are generated by the negative feedback model
90
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
a _ ¼ described as M bM; E_ ¼ cM dE through Bendixson 1 þ KEm criterion. Kurosawa et al. (2002) expanded Griffith's model to the generalized case. They proved that no oscillations are produced by the general class of negative feedback loop described as follows: dx1 ¼ f ðx2 Þ hðx1 Þ; dt
ð1aÞ
dx2 ¼ gðx1 Þ kðx2 Þ; dt
ð1bÞ
where f ðx2 Þ is a monotonically decreasing function, and gðx1 Þ, hðx1 Þ, and kðx2 Þ are monotonically increasing functions. The proof was given by the construction of the Lyapunov function for Eqs. (1). This model corresponds to the case of a common feedback loop of a single gene through transcription and translation. Suppose x1 to be the amount of messenger RNA and x2 to be that df of protein. is negative, given that proteins decrease the rate of dx2 dg is positive, given that gene expression via negative feedback. dx1 dh 4 0 and proteins are produced via the translation of mRNA. dx1 dk 4 0 should hold because of the degradation of the mRNA and dx2 proteins. These conditions should also be satisfied in the simple case of a predator–prey relationship. When x1 and x2 represent the df populations of prey and predator, respectively, o 0 and dx2 dg dh dk 4 0 hold owing to predation, and 4 0 and 4 0 hold dx1 dx1 dx2 because of their death. In the negative feedback with n variables as follows: dx1 ¼ f 1 ðxn Þ h1 ðx1 Þ; dt dx2 ¼ f 2 ðx1 Þ h2 ðx2 Þ; dt
ð2aÞ
ð2bÞ
⋮ dxn ¼ f n ðxn 1 Þ hn ðxn Þ; dt the necessary condition for instability is β 1 β 2 ⋯β n π n 4 sec ; α 1 α 2 ⋯α n n
ð2cÞ
ð3Þ
dhi ðxi Þ df ðx Þ ði ¼ 1; 2; …; nÞ, β i ¼ i i 1 ði ¼ 2; 3; …; nÞ and dxi dxi 1 df ðx Þ β1 ¼ 1 n . This condition is known as the secant condition dxn (Tyson and Othmer, 1978; Thron, 1991; Sontag, 2006). In particular,
where αi ¼
if n ¼2, any values of αi and βi cannot satisfy the condition (3) because the right-hand side of Eq. (3) is infinity. This result is consistent with Kurosawa et al. (2002). The proof that negative feedback with two variables cannot generate oscillation is limited to ordinary differential equations (ODEs). Some researchers have recently focused on cases where ODEs are not applicable. For example, gene expression noise is unavoidable, owing to the low copy numbers of most genes. The effect of molecular noise on the oscillation was assessed via stochastic simulation (Barkai and Leibler, 2000; Gonze et al., 2002; Nishino et al., 2013). Although impossibility of oscillations was proved for ODE model of negative feedback with two variables, whether the other model for negative feedback with
two variables can show self-sustained oscillations has not yet been investigated. Besides ODEs, discrete and ultradiscrete systems are used to model dynamical systems. A discrete system treats the independent variables such as time and space as discrete numbers and is represented as a difference equation (Elaydi, 2005). An ultradiscrete system treats both independent and dependent variables as discrete numbers (Tokihiro et al., 1996). The qualitative behaviors may not generally be conserved by discretizing or ultradiscretizing ODEs. For example, although the logistic equation described as an ODE is attracted to an equilibrium point, the logistic map obtained by the forward discretization of the logistic equation shows a chaotic behavior (May, 1976). In the present study, we focus on discrete and ultradiscrete negative feedback models with two variables. These models can oscillate under some parameter conditions. We also investigated the types of bifurcation and parameter conditions for the oscillations.
2. Discretization In this section, we discretize a negative feedback model with two variables and study its bifurcation structure. 2.1. Model We consider the negative feedback model with two variables, which is reduced from the three-variable model proposed by Kurosawa et al. (2002), as follows: dM ¼ dt
k n aM; P 1þ h
ð4aÞ
dP ¼ sM vP: dt
ð4bÞ
If Eqs. (4) represent a genetic feedback loop, then M and P correspond to amounts of mRNA and protein, respectively. We assume that the parameters k, a, s, v, n, and h are positive. Because Eqs. (4) are included in the class of (1), Eqs. (4) have no limit cycle solutions. To nondimensionalize Eqs. (4), we define the dimenk ks 1 sionless variables M ¼ x, P ¼ 2 y, t ¼ τ and dimensionless v v v 2 a hv parameters b ¼ and c ¼ . Then, Eqs. (4) are transformed into v ks dx 1 ¼ ð5aÞ y n bx; dτ 1þ c dy ¼ x y: dτ
ð5bÞ
By applying tropical discretization (Murata, 2013) to Eqs. (5), we obtain 1 y n bxT 1þ T xT þ 1 xT c ¼ ; δ 1 þ bδ yT þ 1 yT
δ
¼
xT yT ; 1þδ
ð6aÞ ð6bÞ
or, xT þ xT þ 1 ¼
δ
1þ
y n T
c
1 þ bδ
;
ð7aÞ
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
91
Fig. 1. Numerical simulations on the discrete system (Eqs. (7)); (a) n¼2, b¼ 1, c ¼ 0.01 and δ ¼ 3, (b) n¼2, b¼ 1, c ¼0.01 and δ ¼ 1.
yT þ 1 ¼
yT þ δxT ; 1þδ
ð7bÞ
where δ is the time interval. We see that this discretization is relevant because if we substitute τ ¼ T δ and take the continuous limit δ-0, we obtain Eqs. (5). The numerical simulation confirms that the discrete model (7) appears to oscillate or to converge to an equilibrium point depending on parameters n, b, c, and δ (Fig. 1).
From Eq. (11) and b; δ 40, the parameters of n and b, and xn must satisfy n
n 1 4 nbx :
ð12Þ
Note that when δ is sufficiently large, the discrete model (7) has a limit cycle as long as n, b, and xn satisfy Eq. (12). Given that n, b, and xn are positive, parameter n must satisfy n 4 1:
2.2. Stability
We next consider the oscillatory parameter region of b and c.
The parameter region for oscillation can be calculated by stability analysis. Eqs. (7) have a unique fixed point, ðxn ; yn Þ, which satisfies n
bx ¼
ð13Þ
1 n n ; y 1þ c
ð8aÞ
n
Combining Eqs. (8a) and (8b), we obtain bx þ
n
ðbx Þn þ 1 ¼ 1; this is ðbcÞn
rewritten as bc ¼
" #1=n n ðbx Þn þ 1 : n 1 bx
ð14Þ n
xn ¼ yn :
ð8bÞ
Since 0 o
1þ
1yn n o 1, the range of value of xn is c
1 0 o xn o : b
ð9Þ
The Jacobian at ðxn ; yn Þ is 2 3 n 1 nbδ1ð1þbδbx Þ 1 þ bδ 4 5: A¼ δ 1 1þδ
ð10Þ
1þδ
þ bÞδ and The trace and the determinant of A are τ ¼ ð12þþbð1δÞð1 þ δÞ þ nbδ ð1 bx Þ Δ ¼ 1 ð1 , þ bδÞð1 þ δÞ n
2
pffiffiffiffiffiffiffiffiffiffiffiffi
respectively,
and
the
eigenvalues
are
λ ¼ τ 7 τ2 4Δ. 2
If τ2 4Δ o 0, the eigenvalues are complex. τ 2 4Δ is negative for a broad range of values of b and δ, as can be confirmed numerically (Appendix A). In particular, if b is nearly unity, then τ2 4Δ is always negative. We now consider the case of τ2 4Δ o 0. pffiffiffiffi When j λ j ¼ Δ 4 1, i.e., ½n 1 nbx bδ 4 1 þ b; n
From the oscillatory condition of Eq. (12), we obtain bx o n n 1. Combining them, the parameters b and c under oscillatory conditions must satisfy
ð11Þ
the fixed point loses its stability, and the discrete model (7) generates self-sustained oscillations via Neimark–Sacker bifurcation. We can derive that τ2 4Δ o 0 is always satisfied when Eq. (11) holds (Appendix B); i.e., self-sustained oscillation occurs via Neimark–Sacker bifurcation whenever the fixed point is destabilized.
bc o
ðn 1Þðn þ 1Þ=n : n
ð15Þ
Thus, when parameters n, b, c, and δ satisfy Eqs. (11), (13), and (15), the fixed point ðxn ; yn Þ is unstable, and the discrete system (7) behaves as an oscillator. We next consider the range of the value of bc. Eq. (15) can be transformed into bc on1=n ð1 1nÞ1 þ ð1=nÞ . Moreover, n1=n ð1 1nÞ1 þ ð1=nÞ-1 as n-1. Hence, the value of bc is bounded because of the continuity of
ðn 1Þðn þ 1Þ=n . n
ðn 1Þðn þ 1Þ=n is below n 2 hv c ¼ ks , parameters
av h o 1:2: ks
By numerical computation, the maximum of
1.2. Recalling that b and c are defined by b ¼ av and a, v, k, s, and h satisfy ð16Þ
Recall that a and v denote the degradation rates of M and P, respectively, k and s are the synthesis rates of M and P, and h is the threshold for negative regulation. Thus, Eq. (15) suggests that selfsustained oscillation occurs when the degradation rate is lower than the synthesis rate and the threshold of feedback regulation is small. Note that the continuous limit of the Jacobian, limδ-0 A δ I, where I is the identity matrix, equals the Jacobian of the differential equations (5). This result means that if δ is sufficiently small, the period of the discrete system (7) equals that of the continuous system (5).
92
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
ultradiscretization method to Eqs. (22). Let us introduce
2.3. Generalization We will find a general condition for oscillation in discrete negative feedback model with two variables, regardless its function form. We consider the general form of the negative feedback with two variables, that is, Eqs. (1) substituted x1 ¼ x and x2 ¼ y as follows: dx ¼ f ðyÞ hðxÞ; dt
ð17aÞ
dy ¼ gðxÞ kðyÞ; dt
ð17bÞ dk 4 0, dy
where and respectively. By applying tropical discretization to Eqs. (17), we obtain the discrete model xT þ δf ðyT Þ ; xT þ δhðxT Þ
ð18aÞ
y þ δgðxT Þ yT þ 1 ¼ yT T : yT þ δkðxT Þ
ð18bÞ
The Jacobian at the fixed point ðxn ; yn Þ is 2 3 δ xn α 1 δ xn β 1 n 6 1 xn þ δhðxn Þ x þ δhðxn Þ 7 6 7 ; A¼6 n δy β 2 δyn α2 7 4 5 1 n n n n y þ δkðy Þ y þ δkðy Þ
ð19Þ
b¼e
;
yT ¼ eY T =ε ;
c ¼ eC=ε ;
δ ¼ eD=ε :
Then, we obtain
X T þ 1 ¼ X T þ ε log eX T =ε þ eD=ε
δ2 xn yn β1 β2 ; ðxn þ δhðxn ÞÞðyn þ δkðyn ÞÞ
ð24aÞ
X T þ 1 ¼ X T þ max½X T ; max½0; nðY T CÞ þ D max½X T ; X T þ B þ D;
ð25aÞ
Y T þ 1 ¼ Y T þ max½Y T ; X T þ D max½Y T ; Y T þD;
ð25bÞ
where we use the formula for ultradiscretization as follows: lim ε log ½eA1 =ε þeA2 =ε þ ⋯ þ eAn =ε ¼ max½A1 ; A2 ; …; An :
ε- þ 0
ð20Þ
Simulating the ultradiscretized model (25), we can observe oscillations when n Z1, D 4 max½ B; 0, and B þ C o 0 (Fig. 2). These conditions correspond to the ultradiscretization of the oscillatory conditions (11) and (15) for the discrete model, as shown below. From Eqs. (11) and (23), we obtain B=ε
We next focus on the ultradiscretization of our developed discrete model. Eqs. (7) can be rewritten as
yT þ 1 ¼ yT
δ
y n T
yT þ δxT : yT þ δyT
;
ð27Þ
D 4 max½ B; 0: Similarly, from Eqs. (15) and (23), we obtain " # ðn 1Þðn þ 1Þ=n B þ C o ε log : n
ð28Þ
ð29Þ
Taking the limit ε- þ 0, we obtain the ultradiscretized condition B þ C o 0:
ð30Þ
4. Reduction to a Boolean system
3. Ultradiscretization
xT þ 1 ¼ xT
þ BÞ=ε
þ 1 εlog ½n 1:
Assuming that the time interval is large, we can make the discrete model (7) and ultradiscrete model (25) simpler. Taking the limit δ-1 for the discrete model (7), we obtain 1 xT þ 1 ¼ y n ; b 1þ T c yT þ 1 ¼ xT ;
c xT þ δbxT
n
ð21Þ
If the time interval δ is nonzero value, the negative feedback model with two variables (18) might oscillate when the synthesis rates are larger than the degradation rates. By taking continuous limit δ-0, the right-hand side of Eq. (21) diverges to infinity and any values of α1, α2, β1 and β2 cannot satisfy Eq. (21). This means the continuous negative feedback model (17) cannot the generate oscillations. This result corresponds to the secant condition when n ¼2 (Tyson and Othmer, 1978; Thron, 1991; Sontag, 2006).
1þ
ð26Þ
Taking the limit ε- þ 0, we have the ultradiscretized condition
dh dk and α2 ¼ are the degradation rates of x and y, dx dy df dg respectively, β1 ¼ and β 2 ¼ are the synthesis rates of x and y. dy dx If Δ 4 1, the fixed point is unstable. Therefore, the sufficient condition for instability is
β1 β2 xn þ δhðxn Þ yn þ δkðyn Þ 4 þ 1: α1 α2 δα1 xn δα2 yn
ð24bÞ
Taking the limit ε- þ 0, we have the ultradiscrete equations
4 ε log ½e
where α1 ¼
xT þ
ð23Þ
1 ε log ½eX T =ε þ eX T þ B þ D=ε ; 1 þ eðnðY T CÞÞ=ε
D 4 ε log ½e B=ε þ 1 εlog ½n 1 neðX
and the determinant of A is δ xn α 1 δyn α2 1 zw Δ ¼ 1 n x þ δhðxn Þ yn þ δkðyn Þ þ
B=ε
Y T þ 1 ¼ Y T þ ε log ½eY T =ε þ eðX T þ DÞ=ε ε log ½eY T =ε þ eðY T þ DÞ=ε :
df o 0, dh 4 0, dg 40 dy dx dx
xT þ 1 ¼ xT
xT ¼ eX T =ε ;
ð22aÞ
ð22bÞ
Tokihiro et al. (1996) proposed a method for deriving ultradiscrete equations from discrete equations. We will apply the
ð31aÞ
ð31bÞ
where the right-hand side of Eq. (31a) represents a negative feedback effect. According to the simulation of Eqs. (31), this model can display oscillation. Thus, Eqs. (31) represent the simplest discrete model for self-sustained oscillations produced by a negative feedback loop. Similarly, taking the limit D-1 for the ultradiscrete model (25) X T þ 1 ¼ max½0; nðY T CÞ B;
ð32aÞ
Y T þ 1 ¼ XT ;
ð32bÞ
where the right-hand side of Eq. (32a) represents a negative feedback effect. The numerical simulation of Eqs. (32) confirms that this model, as well as the discrete system, conserves the oscillatory solution. The conditions of δ-1 and D-1 can be interpreted as a sufficiently long time interval between successive reactions in
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
93
Fig. 2. Numerical simulations on the ultradiscrete system (Eqs. (25)); (a) n¼ 1, B ¼0, C ¼ 1 and D ¼1, (b) n ¼1, B¼ 1, C ¼ 1 and D¼ 1.
biochemical systems or as cycles between long generation times in ecosystems. These simple models generally correspond to negative feedback loops where the reaction interval is relatively larger than reaction rate. Moreover, Eqs. (32) can be reduced to a Boolean system under the condition shown below. A Boolean system is one that consists of binary variables and logical operations, such as AND ( 4 ), OR ( 3 ), and NOT (:). Suppose n ¼1, C r 0, and C B ¼ 1. Eqs. (32) become X T ¼ :Y T ;
ð33aÞ
Y T ¼ XT ;
ð33bÞ
where 1 Y T is replaced with :Y T . This model is a Boolean system because the range of XT and YT is restricted in f0; 1g if we assign 0 or 1 to the variables as initial values. This means that the ultradiscrete model for negative feedback loops, Eqs. (25), includes a Boolean system, Eqs. (33), as a special case.
5. Discussion and conclusions We focused on the simple negative feedback model that consists of two elements. As Tyson and Othmer (1978), Thron (1991) and Kurosawa et al. (2002) proved, the ODE model with two variables never oscillates, regardless of its function form. We found that the discretized and ultradiscretized models can oscillate depending on the values of the parameters. In discrete model, we derived the general condition for oscillation, which covers any forms of function. Discretizing and ultradiscretizing can contribute to destabilizing a stable negative feedback system. We have also found the developed ultradiscrete model can be reduced to a Boolean system under some conditions. Some textbooks illustrate the mechanism of circadian rhythms produced by a simple negative feedback loop as follows: (i) expression of mRNA induces the production of proteins, (ii) proteins inhibit the expression via feedback, (iii) lowered expression level decreases protein amounts, (iv) the expression starts again with the disappearance of the inhibitory proteins (Alberts et al., 2008). This model is equivalent to the Boolean system Eqs. (33). Note that a discretized or ultradiscretized simple feedback model can show oscillations but that a continuous limit of the model cannot always conserve its rhythmicity. This observation may explain why most biological rhythms consist of many molecular parts. The parameter n represents the nonlinearity of negative feedback control. n corresponds to an index of the cooperativity of protein–DNA interaction in the context of a gene regulatory network. Some studies have shown that higher cooperativity is more likely to contribute in the destabilization of the negative
feedback system described in an ODE (Leloup and Goldbeter, 1998; Kurosawa et al., 2002). For example, the negative feedback model with three variables proposed in Kurosawa et al. (2002) can be destabilized when n Z 9. In contrast, our developed model can be destabilized when n 4 1 for the discrete model and n Z 1 for the ultradiscrete model, suggesting that the nonlinearity of feedback is not necessary for discretized oscillators. We used the tropical discretization for discretizing Eqs. (5). This discretization was proposed by Murata (2013) as a systematic manner of discretizing and ultradiscretizing first-order differential equations. This discretization preserves the positivity of variables if the initial values and parameters are positive. Moreover, dependent variables x and y are proportional to the amounts of some proteins or the populations of prey and predator. Thus, x and y must be positive, and tropical discretization is appropriate in this case. In contrast, applying the other discretization to Eqs. (5), the system may not preserve positivity. For example, applying the forward difference to Eqs. (5), the discrete model has a negative value when δ is large. We consider the biological interpretation of the tropical discretization. The right-hand sides of Eqs. (6) are the right-hand sides of Eqs. (5) divided by 1 þ bδ or 1 þ δ. Thus, the larger the time-step size δ, the smaller are the right-hand sides of Eqs. (6). This result means that when the time interval between a chemical reaction and the next chemical reaction is large, the effect of synthesis and the degradation of a product are small. This assumption may be valid for real biological systems to preserve the positivity of components. As mentioned in the Introduction, the stochastic model can be an alternative to reproduce the behavior of dynamical systems in nature. We transformed the deterministic model (4) into a stochastic model (Appendix C). The obtained model produced pulses intermittently for a broad range of the values of parameters. This occurs because the state drifting around the fixed point happens to be ejected due to the molecular noise and randomly makes a large-amplitude rotation. Unlike the discrete and ultradiscrete model, we did not succeed in observing any oscillations with a constant period and amplitude. The discretized or ultradiscretized model corresponds to some actual biological systems. Behavior as a discretized system is often essential to an ecological system. Many species of plants and animals perform reproductive activity every year. The population dynamics can be described as a discrete system with δ ¼1 [year]. It is the common case that the production rate of species Y is enhanced by species X and conversely the species Y decrease the production rate of species X. Based on our research, we can conclude that this simple ecological negative feedback regulation can cause the oscillation as long as the condition (21) is satisfied. The number of population should be represented as a discrete value rather than continuous one when the distribution of individuals is sparse and
94
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
Appendix A. Condition for τ2 4Δ o 0
12
We numerically computed the sign of τ2 4Δ, which is the discriminant of the characteristic polynomial of the Jacobian shown in Eq. (10) (Fig. 3). We found that τ2 4Δ is negative for most values of b; δ. In particular, τ2 4Δ is always negative when b is near b¼1, resulting in self-sustained oscillation or damped oscillation. Recall that b is defined as b ¼ av and that a and v are degradation rates of M and P, respectively. Hence, b ¼1 means that the degradation rate of M nearly equals that of P.
10 8 6 4
Appendix B. The condition for τ2 4Δ o 0 includes the condition for Δ 41
2 0
0
2
4
6
8
10
12
Numerical calculation suggests that Δ 4 1 is satisfied when τ2 4Δ o 0 holds and when n ¼ 3; c ¼ 0:8 (Fig. 3). Moreover, the area of τ2 4Δ o0 includes the area of Δ 4 1 for arbitrary n and δ. This is analytically confirmed as follows:
Number of
Fig. 3. Sign of τ2 4Δ. The black area represents the sets of b and c that satisfy τ2 4Δ 40 when n¼ 3 and c¼ 0.8. The broken line denotes Δ ¼ 1 at which Neimark– Sacker bifurcation occurs.
τ2 4Δ ¼ ½ð2 þ ð1 þ bÞδÞ2 4f1 þ nbδ2 ð1 bxn Þgð1 þ bδÞð1 þ δÞ =ð1 þbδÞ2 ð1 þ δÞ2
200
¼ ½ 4nb ð1 bx Þδ 4nbð1 þbÞð1 bx Þδ
100
Assuming Δ 41, we obtain δ 4
0
n
2
n
2
4nbð1 bx Þ þ ð1 bÞ δ =ð1 þ bδÞ ð1 þ δÞ2 : n
0
20
40
60
Time Fig. 4. Stochastic simulations on the system Eqs. (25); h ¼ a ¼ v ¼ 1, k ¼ v ¼ 5, Ω¼ 100 and τ¼ 1.
ðB:1Þ
2
2
2
ðB:2Þ
1þb . By substituting the n ½n 1 nbx b inequality into the second term of the numerator in Eq. (B.2), we have 1þb τ2 4Δ o 4nb2 ð1 bxn Þδ2 4nbð1 þ bÞð1 bxn Þ n ½n 1 nbx b i 2 n 4nbð1 bx Þ þ ð1 bÞ2 δ =ð1 þ bδÞ2 ð1 þ δÞ2 ðB:3Þ "
ð1 bÞ2 n n 1 nbx n
nð3b þ 1Þðb þ 3Þð1 bx Þ 2 δ =ð1 þ bδÞ2 ð1 þ δÞ2 : þ n n 1 nbx
¼ 4nbð1 bx Þðbδ þ 1Þ þ n
the number of species is quite low. Ultradiscrete method should be applicable in this case. At the molecular level, biochemical reaction can be also described using the discrete model due to time delay. For example, intron produces the time delay in a transcriptional process. It is theoretically well known that the time delay can contribute to destabilize the system. Some experimental research confirms this hypothesis. Hes7 gene that is responsible for segmentation clock forms autonegative feedback loop, and shows the self-sustained oscillations of gene expression. The intronic region of Hes7 is required for the genetic oscillations (Takashima et al., 2011; Harima et al., 2013). In addition, the length of intron can lengthen the period of transcriptional oscillation in engineered genetic oscillator (Swinburne et al., 2008). Meanwhile, post-translational modification and nuclear transportation of proteins can generate time delay in circadian rhythms. The discrete model can be used to explicitly incorporate the time delay between synthesis and being mature molecule that can function. The ultradiscrete method should be used when the number of molecules is low as well as the ecological level. The effect of a low number of individuals or molecules needs to be addressed through an ultradiscrete model in future studies.
n
2
ðB:4Þ
n
Because 1 bx 4 0 and n 1 nbx 4 0 can be derived from Eqs. (9) and (12), respectively, we conclude τ2 4Δ o 0 if Δ 4 1. Appendix C. Stochastic simulation for genetic feedback loop with two variables We checked the rhythmicity of the stochastic model transformed from the deterministic feedback model with two variables. If Eqs. (4) represent a genetic feedback loop, the ODE model can be decomposed into some elementary biochemical reactions that consist of a gene G, mRNA M and protein P as follows: G þ nP-GP n ;
ðC:1aÞ
GP n -G þ nP;
ðC:1bÞ
G-G þM;
ðC:1cÞ
M-∅;
ðC:1dÞ
M-M þ P;
ðC:1eÞ
Acknowledgment
P-∅:
ðC:1f Þ
We would like to thank T. Kawabe (Kyushu University) for many helpful discussions. This research was supported by Inamori Foundation and Takeda Science Foundation (H.I.).
The two tional to time of (Nishino
scale parameters Ω and τ are introduced; Ω is proporthe number of molecules in the system, and τ scales the association and dissociation of proteins to the DNA et al., 2013). We assumed the transition rates of the
S. Gibo, H. Ito / Journal of Theoretical Biology 378 (2015) 89–95
reactions (C.1a)–(C.1f) are 1τ ðPhΩÞn G, 1τ ð1 GÞ, kΩG, aM, sM and vP. We performed stochastic simulation for the reactions using the Gillespie algorithm and observed stochastic pulse for a broad range of parameters (Fig. 4).
References Alberts, B., Johnson, A., Lewis, J., Raff, M., Keith, R., Peter, W., 2008. Molecular Biology of the Cell, 5th ed. Garland Science, New York. Alon, U., 2007. Network motifs: theory and experimental approaches. Nat. Rev. Genet. 8 (6), 450–461. Aronson, B.D., Johnson, K.A., Loros, J.J., Dunlap, J.C., 1994. Negative feedback defining a circadian clock: autoregulation of the clock gene frequency. Science 263, 1578–1584. Barkai, N., Leibler, S., 2000. Biological rhythms: circadian clocks limited by noise. Nature 403, 267–268. Elaydi, S., 2005. An Introduction to Difference Equations, 3rd ed. Springer, New York. Fussmann, G.F., Ellner, S.P., Shertzer, K.W., Hairston, N.G., 2000. Crossing the Hopf bifurcation in a live predator–prey system. Science 290, 1358–1360. Gonze, D., Halloy, J., Goldbeter, A., 2002. Robustness of circadian rhythms with respect to molecular noise. Proc. Natl. Acad. Sci. USA 99, 673–678. Griffith, J.S., 1968. Mathematics of cellular control processes I. Negative feedback to one gene. J. Theor. Biol. 20, 202–208. Harima, Y., Takashima, Y., Ueda, Y., Ohtsuka, T., Kageyama, R., 2013. Accelerating the tempo of the segmentation clock by reducing the number of introns in the Hes7 gene. Cell Rep. 3, 1–7. Hirata, H., Yoshiura, S., Ohtsuka, T., Bessho, Y., Harada, T., Yoshikawa, K., Kageyama, R., 2002. Oscillatory expression of the bHLH factor Hes1 regulated by a negative feedback loop. Science 298, 840–843. Ishiura, M., Katsuna, S., Aoki, S., Iwasaki, H., Andersson, C.R., Tanabe, A., Golden, S.S., Johnson, C.H., Kondo, T., 1998. Expression of a gene cluster kaiABC as a circadian feedback process in cyanobacteria. Science 281, 1519–1523. Kurosawa, G., Mochizuki, A., Iwasa, Y., 2002. Comparative study of circadian clock models, in search of processes promoting oscillation. J. Theor. Biol. 216, 193–208.
95
Lahav, G., Rosenfeld, N., Sigal, A., Geva-Zatorsky, N., Levine, A.J., Elowitz, M.B., Alon, U., 2004. Dynamics of the p53-Mdm2 feedback loop in individual cells. Nat. Genet. 36, 147–150. Leloup, J.C., Goldbeter, A., 1998. Modeling circadian oscillations of PER and TIM proteins in Drosophila. In: Touitou, Y. (Ed.), Biological Clocks. Mechanisms and Applications. Elsevier Science, Amsterdam, pp. 81–88. Liebhold, A., Elkinton, J., Williams, D., Muzika, R.M., 2000. What causes outbreaks of the gypsy moth in North America? Popul. Ecol. 42, 257–266. May, R.M., 1976. Simple mathematical models with very complicated dynamics. Nature 261, 459–467. McWatters, H.G., Devlin, P.F., 2011. Timing in plants—a rhythmic arrangement. FEBS Lett. 585 (10), 1474–1484. Murata, M., 2013. Tropical discretization: ultradiscrete Fisher-KPP equation and ultradiscrete Allen–Cahn equation. J. Differ. Equ. Appl. 19 (6), 1008–1021. Nishino, R., Sakaue, T., Nakanishi, H., 2013. Transcription fluctuation effects on biochemical oscillations. PLoS ONE 8 (4), e60938. Partch, C.L., Green, C.B., Takahashi, J.S., 2014. Molecular architecture of the mammalian circadian clock. Trends Cell Biol. 24 (2), 90–99. Rosenfeld, N., Elowitz, M.B., Alon, U., 2002. Negative autoregulation speeds the response times of transcription networks. J. Mol. Biol. 323 (5), 785–793. Shen-Orr, S.S., Milo, R., Mangan, S., Alon, U., 2002. Network motifs in the transcriptional regulation network of Escherichia coli. Nat. Genet. 31, 64–68. Sontag, E.D., 2006. Passivity gains and the “secant condition” for stability. Syst. Control Lett. 55, 177–183. Swinburne, I.A., Miguez, D.G., Landgraf, D., Silver, P.A., 2008. Intron length increases oscillatory periods of gene expression in animal cells. Genes Dev. 22, 2342–2346. Takashima, Y., Ohtsuka, T., González, A., Miyachi, H., Kageyama, R., 2011. Intronic delay is essential for oscillatory expression in the segmentation clock. Proc. Natl. Acad. Sci. USA 108 (8), 3300–3305. Thron, C.D., 1991. The secant condition for instability in biochemical feedback control—I. The role of cooperativity and saturability. Bull. Math. Biol. 53 (3), 383–401. Tokihiro, T., Takahashi, D., Matsukidaira, J., Satsuma, J., 1996. From soliton equations to integrable cellular automata through a limiting procedure. Phys. Rev. Lett. 76 (18), 3247–3250. Tyson, J.J., Othmer, H.G., 1978. The dynamics of feedback control circuits in biochemical pathways. Prog. Theor. Biol. 5, 1–63.