Biol Cybern (2014) 108:365–380 DOI 10.1007/s00422-014-0605-7

ORIGINAL PAPER

Fixed-point bifurcation analysis in biological models using interval polynomials theory Gerasimos G. Rigatos

Received: 27 August 2013 / Accepted: 15 April 2014 / Published online: 10 May 2014 © Springer-Verlag Berlin Heidelberg 2014

Abstract The paper proposes a systematic method for fixed-point bifurcation analysis in circadian cells and similar biological models using interval polynomials theory. The stages for performing fixed-point bifurcation analysis in such biological systems comprise (i) the computation of fixed points as functions of the bifurcation parameter and (ii) the evaluation of the type of stability for each fixed point through the computation of the eigenvalues of the Jacobian matrix that is associated with the system’s nonlinear dynamics model. Stage (ii) requires the computation of the roots of the characteristic polynomial of the Jacobian matrix. This problem is nontrivial since the coefficients of the characteristic polynomial are functions of the bifurcation parameter and the latter varies within intervals. To obtain a clear view about the values of the roots of the characteristic polynomial and about the stability features they provide to the system, the use of interval polynomials theory and particularly of Kharitonov’s stability theorem is proposed. In this approach, the study of the stability of a characteristic polynomial with coefficients that vary in intervals is equivalent to the study of the stability of four polynomials with crisp coefficients computed from the boundaries of the aforementioned intervals. The efficiency of the proposed approach for the analysis of fixed-point bifurcations in nonlinear models of biological neurons is tested through numerical and simulation experiments. Keywords Fixed-point bifurcation analysis · Nonlinear biological neurons · Characteristic polynomials · Interval polynomials · Kharitonov’s stability theory G.G. Rigatos (B) Unit of Industrial Automation, Industrial Systems Institute, Stadiou str, 26504Rion Patras, Greece e-mail: [email protected]

1 Introduction The paper proposes the use of interval polynomial theory for the study of fixed-point bifurcations and the associated stability characteristics in circadian cells and similar biological models. The bifurcation parameter is usually one of the coefficients of the biological system or a feedback control gain that is used to modify the system’s dynamics. The bifurcation parameter varies within intervals; therefore, assessment of the stability features for all fixed points of the system (each fixed point depends on a different value of the bifurcation parameter) requires extended computations. Therefore, it is significant to develop methods that enable to draw a conclusion about the stability features of larger sets of fixed points where each set is associated with an interval from which the bifurcation parameter takes its values. A solution to this problem can be obtained from interval polynomials theory and Kharitonov’s theorem (Rigatos and Siano 2011; Toscano and Lyonnet 2010). Aiming at understanding how the dynamics of neurons is modified due to external stimuli and parametric variations, the topic of fixed-point bifurcations in such biological models has been widely studied (Lysetskiy and Zurada 2004; Wang et al. 2013; Xu et al. 2006]. Similar results and methods have been extended to artificial neural networks, such as Hopfield associative memories (Haschke and Steil 2005; Song and Wei 2005). Additionally, studies on bifurcations appearing in biological models such as circadian cells show how the dynamics and the stability of such models can be affected by variations of external inputs and internal parameters (Nagy 2008; Tsumoto et al. 2006). A subject of particular interest has been also the control of bifurcations appearing in biological systems which is implemented with state feedback and which finally succeeds to modify the biological oscillator’s dynamics and its convergence to attractors (Dorak

123

366

2013; Duan et al. 2012; Efimov and Fradkov 2006; Vidal et al. 2012; Nguyen and Hong 2012; Nguyen et al. 2012; Wang et al. 2004). In this paper, as a case study, bifurcations of fixed points in two particular types of biological models will be examined. The first model is the FitzHugh–Nagumo neuron. The voltage of the neuron’s membrane exhibits oscillatory variations after receiving suitable external excitation either when the neuron is independent from neighboring neural cells or when the neuron is coupled to the neighboring neural cells through synapses or gap junctions. The FitzHugh–Nagumo model of biological neurons is a simplified model (2nd order) of the Hodgkin–Huxley model (4th order), where the latter is considered to reproduce efficiently in several cases the neuron’s dynamics. The FitzHugh–Nagumo model describes the variations of the voltage of the neuron’s membrane as a function of the ionic currents that get through the membrane and of an external current that is applied as input to the neuron (Adil et al. 2012; Baird and Terman 2009; Rigatos and Rigatou 2013; Zhang et al. 2007). The second model is that of Neurospora circadian cells. Circadian cells perform proteins synthesis within a feedback loop. The concentration of the produced proteins exhibits periodic variations in time. The circadian cells regulate the levels of specific proteins’ concentration through a RNA transcription and translation procedure (Laroche and Claude 2004; Leloup et al. 2006; That and Ding 2012). Typically, the stages for analyzing the stability features of the fixed points that lie on bifurcation branches comprise (i) the computation of fixed points as functions of the bifurcation parameter and (ii) the evaluation of the type of stability for each fixed point through the computation of the eigenvalues of the Jacobian matrix that is associated with the system’s nonlinear dynamics model . Stage (ii) requires the computation of the roots of the characteristic polynomial of the Jacobian matrix. This problem is nontrivial since the coefficients of the characteristic polynomial are functions of the bifurcation parameter and the latter varies within intervals (Chen et al. 2000; Liao et al. 2001; Nagy 2009). To obtain a clear view about the values of the roots of the characteristic polynomial and about the stability features they provide to the system, the use of interval polynomials theory and particularly of Kharitonov’s stability theorem is proposed (Rigatos and Siano 2011; Toscano and Lyonnet 2010). In this approach, the study of the stability of a characteristic polynomial with coefficients that vary in intervals is equivalent to the study of the stability of four polynomials with crisp coefficients computed from the boundaries of the aforementioned intervals. The significance of the proposed stability analysis method for biological systems is easy to be understood. Processes and phenomena in biological systems either (i) exhibit asymptotic stability properties and convergence to stable fixed points which means that the effects of short-term excitation

123

Biol Cybern (2014) 108:365–380

become extinct or that the transient effects of long-term external inputs progressively attenuate and at the same time specific system parameters reach levels that cannot be further changed or (ii) exhibit limit cycles which means that periodic and sustained-in-time variations appear. These oscillatory phenomena are frequently met in hormone and protein synthesis cycles, in the hematopoietic process, in the functioning of body’s organs, etc. Therefore, the biology of what is done in the paper has to do with tools that enable to assess the status of biological systems in terms of stability properties and consequently in terms of closeness to healthy conditions. The loss of stability and periodicity in biological systems is associated with the appearance of diseases. By identifying which are the parameters or feedback gains in the dynamic model of biological systems which are responsible for the loss of stability or periodicity, it is possible to decode diseases’ mechanisms and to improve the effects of pharmaceutical schemes. The structure of the paper is as follows: In Sect. 2, the biological significance of FitzHugh–Nagumo neurons and of circadian oscillators is explained. In Sect. 3, the generalization of Routh–Hurwitz criterion through Kharitonov’s theorem for dynamical systems subjected to parametric uncertainty is explained. In Sect. 4, the stages in analyzing bifurcations in dynamical systems are presented. In Sect. 5, bifurcation analysis for the FitzHugh–Nagumo model of biological neurons is given. In Sect. 3, bifurcation analysis in circadian cells is given. In Sect. 7, the problem of feedback control of bifurcations is analyzed, and it is shown that the feedback control gains can function as bifurcating parameters that determine the type of stability of the system’s fixed points. In Sect. 8, simulation tests are provided to further confirm the results of the previous fixed-point bifurcation and stability analysis. Finally, in Sect. 9, concluding remarks are stated.

2 The FitzHugh–Nagumo and the circadian oscillators’ model 2.1 The FitzHugh–Nagumo model of neuron oscillations The FitzHugh–Nagumo neuron is a simplified pointwise model of voltage variations in the neuron’s membrane. Voltage dynamics is best described by the cable’s partial differential equation or wave-type partial differential equations, and this model provides the spatiotemporal voltage variations along the membrane (Baird and Terman 2009; Gerstner and Kistler 2002; Rigatos 2013). To avoid the use of this computationally demanding model, pointwise study of voltage variations with respect to time has been considered and this has resulted in the Hogdkin–Huxley model. This is a model of four coupled nonlinear differential equations.

Biol Cybern (2014) 108:365–380

367

To further reduce the complexity of the Hodgkin–Huxley model, the FitzHugh–Nagumo model of the neuron’s dynamics has been introduced, which comprises only two coupled differential equations. The FitzHugh–Nagumo model is a two-dimensional simplified form of the Hogdkin–Huxley model and in several cases represents efficiently the membrane’s voltage dynamics in neurons. The model captures all the qualitative properties of the Hodgkin–Huxley dynamics, in the sense that it results in an equivalent bifurcation diagram (loci of fixed points with respect to model’s parameters).

3 Generalization of the Routh–Hurwitz criterion with Kharitonov’s theorem

2.2 The model of circadian oscillators

p(λ) = cn λn + cn−1 λn−1 + · · · + c1 λ + c0

Circadian neurons are a small cluster of 2 × 104 neurons in the ventral hypothalamus and provide the body with key time-keeping signals, known as circadian rhythms. Circadian neurons mainly receive input from the retina, can be coupled through neurotransmitters and can stimulate other cells in the body also exhibiting an oscillatory behavior that provides time keeping for other functions of the body. There are also peripheral circadian oscillators which are responsible for rhythmicity in the functioning of various organs in the human body. The peripheral circadian oscillators are controlled by the circadian neurons through messages sent by the latter, in the form of hormones (ACTH, cortisol, etc.), metabolites, inputs from the sympathetic neural system or variations in body temperature (François 2005; Gonze et al. 2000, 2003; Leloup et al. 2006). Circadian rhythms affect cell cycle apoptosis and DNA repair, drug metabolism and detoxification, as well as angiogenesis. The disruption of periodicity and synchronism in circadian oscillators is known to be responsible for the appearance of malignancies (Clairambault 2006, 2008; Lévi 2008). On the other hand, pharmacological treatment of tumors when taking into account the variations of proteins concentration controlled by the circadian oscillators can be more efficient. Mathematical models of circadian oscillators consist of sets of ordinary differential equations describing the feedback loop of the generation of proteins in the circadian cells. The main stages of these feedback loops are transcription of genes into mRNA, translation of mRNA into proteins and inhibition of the transcription stage by the generated protein. The number of ordinary differential equations in such models is moderate in simple organisms (e.g., neurospora and drosophila), whereas it becomes increased in the case of mammals. In this paper, a Kharitonov theory-based generalization of the Routh–Hurwitz stability criterion will be used for bifurcation analysis and assessment of the stability properties in the aforementioned models of biological systems. The basic stages of this approach are explained next.

The Routh–Hurwitz criterion gives a systematic method for finding if all the roots of a polynomial lie in the left complex semi-plane and thus are stable. Moreover, one can find out the number of complex eigenvalues which lie in the right complex semi-plane and thus are unstable. Additionally, one can find if there are roots of the polynomial with purely imaginary part. The Routh–Hurwitz stability analysis method is outlined in the “Appendix.” Next, the following characteristic polynomial is considered (1)

The polynomial’s coefficients vary within intervals, that is cimin ≤ci ≤ci max i = 0, 1, 2, · · · , n

(2)

The study of the stability of the above characteristic polynomial is now transferred to the study of the stability of four Kharitonov polynomials, where each polynomial is written as i λn−1 +· · ·+c1i λ+c0i i = 1, · · · , 4 (3) pi (λ) = cni λn +cn−1

The associated coefficients are written as i i , c2k+1 k = m, m − 1, · · · , 0, 1 c2k

(4)

where m = n/2 if n is an even number and m = (n − 1)/2 if n is an odd number. For the first Kharitonov polynomial p1 (s) = cn1 s n + 1 cn−1 s n−1 + · · · + c11 s + c01 , the coefficients are 1 max = c2k if k is an even number c2k 1 max c2k+1 = c2k+1 if k is an even number 1 min c2k = c2k if k is an odd number

(5)

1 min c2k+1 = c2k+1 if k is an odd number

For the second Kharitonov polynomial p2 (s) = cn2 s n + + · · · + c12 s + c02 , the coefficients are

2 s n−1 cn−1

2 min c2k = c2k if k is an even number 2 min c2k+1 = c2k+1 if k is an even number 2 max c2k = c2k if k is an odd number

(6)

2 max c2k+1 = c2k+1 if k is an odd number

123

368

Biol Cybern (2014) 108:365–380

For the third Kharitonov polynomial p3 (s) = cn3 s n + + · · · + c13 s + c03 , the coefficients are

3 s n−1 cn−1

3 min = c2k c2k 3 c2k+1 3 c2k 3 c2k+1

= = =

max c2k+1 max c2k min c2k+1

if k is an even number if k is an even number if k is an odd number

(7)

if k is an odd number

For the fourth Kharitonov polynomial p4 (s) = cn4 s n + + · · · + c14 s + c04 , the coefficients are

4 s n−1 cn−1

4 max = c2k c2k 4 c2k+1 4 c2k 4 c2k+1

= = =

min c2k+1 min c2k max c2k+1

if k is an even number if k is an even number if k is an odd number

(8)

if k is an odd number

As an example, the following characteristic polynomial is considered p(λ) = λ3 +(a+b)λ2 +abλ+ K with parametric uncertainty given by K ∈[8, 12], a∈[2.5, 3.5], b∈[3.8, 4.2]. The coefficients of the characteristic polynomial have the following variation ranges: c3 = 1, c2 = (a + b)∈[6.3, 7.7], c1 = ab∈[9.5, 14.7] and c0 = K ∈[8, 12]. The associated Kharitonov polynomials are p1 (λ) = 12 + 14.7λ + 6.3λ2 + λ3 p2 (λ) = 8 + 14.7λ + 7.7λ2 + λ3 p3 (λ) = 8 + 9.5λ + 7.7λ2 + λ3 p4 (λ) = 12 + 9.5λ + 6.3λ + λ 2

(9)

(10)

(11)

Since there is no change of sign in the first column of the Routh matrix, it can be concluded that polynomial p2 (λ) has stable roots. For polynomial p3 (λ), λ3 | 1 9.5 λ2 | 7.7 8 λ1 | 8.4610 0 λ0 | 8

123

(13)

Since there is no change of sign in the first column of the Routh matrix, it can be concluded that polynomial p4 (λ) has stable roots. Using that all four Kharitonov characteristic polynomials pi (λ) have stable roots, it can be concluded that the initial characteristic polynomial p(λ) has stable roots when its parameters vary in the previously defined ranges. Equivalently, the computation of the four Kharitonov characteristic polynomials can be performed as follows: p1 (λ) = c0max + c1max λ + c2min λ2 + c3min λ3 + c4max λ4 + c5max λ5 + c6min λ6 + · · · p2 (λ) = c0min + c1min λ + c2max λ2 + c3max λ3 + c4min λ4 + c5min λ5 + c6max λ6 + · · · p3 (λ) = c0min + c1max λ + c2max λ2 + c3min λ3

(14)

p4 (λ) = c0max + c1min λ + c2min λ2 + c3max λ3 + c4max λ4 + c5min λ5 + c6min λ6 + · · ·

3

Since there is no change of sign in the first column of the Routh matrix, it can be concluded that polynomial p1 (λ) has stable roots. For polynomial p2 (λ), 1 14.7 λ3 | 8 λ2 | 7.7 λ1 | 13.6610 0 8 λ0 |

λ3 | 1 9.5 λ2 | 6.3 12 λ1 | 7.5952 0 λ0 | 12

5 max 6 + c4min λ4 + λmax 5 λ + c6 λ + · · ·

Next, Routh–Hurwitz criterion is applied for each one of the four Kharitonov polynomials. For polynomial p1 (λ), λ3 | 1 14.7 λ2 | 6.3 12 λ1 | 12.795 0 λ0 | 12

Since there is no change of sign in the first column of the Routh matrix, it can be concluded that polynomial p3 (λ) has stable roots. For polynomial p4 (λ),

(12)

while an equivalent statement of Kharitonov’s theorem is as follows (Rigatos and Siano 2011; Toscano and Lyonnet 2010): Theorem Each characteristic polynomial in the interval polynomials family p(λ) is stable if and only if all four Kharitonov polynomials pi (λ), i = 1, · · · , 4 are stable.

4 Stages in bifurcations analysis Some terminology associated with fixed points is as follows: Consider the autonomous dynamical system x˙ = Ax. A fixed point for this system is called hyperbolic if none of the eigenvalues of matrix A has zero real part. A hyperbolic fixed point is called a saddle if some of the eigenvalues of matrix A have real parts greater than zero and the rest of the eigenvalues have real parts less than zero. If all of the eigenvalues have negative real parts then the hyperbolic fixed point is called a stable node or sink. If all of the eigenvalues have positive real parts then the hyperbolic fixed point is called an unstable node or source. If the eigenvalues are purely imaginary then one has an elliptic fixed point which is said to be a center.

Biol Cybern (2014) 108:365–380

369

Bifurcation of the equilibrium point means that the locus of the equilibrium on the plane (having as horizontal axis the bifurcation parameter and as vertical axis the fixed-point variable) changes due to variation of the system’s parameters. Equilibrium x ∗ is a hyperbolic equilibrium point if the real parts of all eigenvalues of the Jacobian matrix are nonzero. One can have stable or unstable bifurcations of hyperbolic equilibria. A saddle-node bifurcation occurs when one eigenvalue is zero while the rest of the eigenvalues have nonzero real part. A Hopf bifurcation appears when the hyperbolicity of the equilibrium point is lost due to variation of the system parameters and the eigenvalues become purely imaginary. By changing the values of the parameters at a Hopf bifurcation, an oscillatory solution appears. The stages for finding bifurcations in nonlinear dynamical systems are described next. The following autonomous differential equation is considered:

dx = f 1 (x, q) dt

(15)

where x is the state vector x∈R n and q∈R m is the vector of the system parameters. In Eq. (15), a point x ∗ satisfying f 1 (x ∗ ) = 0 is an equilibrium point. Therefore, from the condition f 1 (x ∗ ) = 0, one obtains a set of equations which provide the equilibrium point as function of the bifurcating parameter. The stability of the equilibrium point can be evaluated by linearizing the system’s dynamic model around the equilibrium point and by computing eigenvalues of the Jacobian matrix. The Jacobian matrix at the equilibrium point can be written as

5 Bifurcation analysis of the FitzHugh–Nagumo neuron 5.1 Equilibria and stability of the FitzHugh–Nagumo neuron The model of the dynamics of the FitzHugh–Nagumo neuron is considered dV = V (V − a)(1 − V ) − w + I dt dw = (V − γ w) dt

(18)

By defining the state variables x1 = w, x2 = V and setting current I as the external control input, one has d x1 = −γ x1 + x2 dt d x2 = x2 (x2 − a)(1 − x2 ) − x1 + u dt y = h(x) = x1

(19)

Equivalently one has       x˙1 −γ x1 + x2 0 = + u x˙2 x2 (x2 − a)(1 − x2 ) − x1 1 y = h(x) = x1

(20)

Next, the fixed points of the dynamic model of the FitzHugh– Nagumo neuron are computed. It holds that x˙1 = −γ x1 + x2 x˙2 = x2 (x2 − α)(1 − x2 ) − x1 + I

(21)

Setting I = 0, and setting x˙1 = 0 and x˙2 = 0 (nullclines), the associated fixed points are found x2 = −γ x1

D f 1 (x ∗ ) =

∂ f 1 (x) |x=x ∗ ∂x

(16)

and the determinant of the Jacobian matrix provides the characteristic equation which is given by

(22)

In this example, the numerical values of the parameters are taken to be γ = 0.2 and a = 7. Then, by substituting the relation for x2 from the first row of Eq. (22) into the second row, one obtains x1 [γ 3 x12 + γ 2 (1 + a)x1 + (aγ − 1)] = 0

det(λi In − D f 1 (x ∗ )) = cn λin + cn−1 λin−1 + · · · + c1 λi + c0 = 0,

x1 (−γ 3 x12 + (aγ 2 + γ 2 )x1 − αγ − 1) = 0

(17)

where In is the n×n identity matrix, and μi with i = 1, 2, · · · , n denotes the eigenvalues of the Jacobian matrix D f 1 (x ∗ ). Using the above characteristic polynomial, one can formulate conditions for the system to have stable or unstable fixed points or saddle-node fixed points. From the requirement, the eigenvalues of the system to be purely imaginary one obtains a condition, i.e., values that the bifurcating parameter should take, so as Hopf bifurcations to appear.

(23)

From the above conditions, one has the first fixed point x1∗,1 = 0 x2∗,1 = 0

(24)

whereas, since the determinant  = γ 4 (1 + a) − 2aγ 4 − 4γ 3 ) < 0, no other real valued solutions for (x1 , x2 ) can be found and consequently no other fixed points can be considered. Next, the type of stability that is associated with the fixed point x1∗,1 = 0, x2∗,1 = 0 is examined. The Jacobian matrices of the right-hand side of Eq. (21) are computed, and the

123

370

Biol Cybern (2014) 108:365–380

associated characteristic polynomials are found     ∂ f1 ∂ f1 −γ  ∂ x ∂ x 1 2 J = ∂ f2 ∂ f2 = −1 −3x22 + 2(1 + a)x2 − a ∂x ∂x 1

(25)

2

Stability at the fixed point x1 = 0 and x2 = 0: Substituting the values of x1 and x2 in Eq. (27), the obtained characteristic polynomial is p(λ) = λ2 + (a + γ )λ + (αγ + ). Next, it is considered that a∈[amin , amax ] = [6.8, 7.2], whereas  = 1.5 and γ = 0.2 are crisp numerical values. Then, for the coefficients of the characteristic polynomial p(s) = λ2 + k1 λ + k2 , one has k1 = (a + γ ) ∈ [7.1, 7.5] and k2 = (aγ ) +  ∈ [3.54, 3.66]. The application of the Routh–Hurwitz criterion gives λ2 | 1 k 2 λ1 | k 1 0 λ0 | k 2

(26)

Since k1 > 0 and k2 > 0, there is no change of sign in the first column of the Routh matrix and the roots of the characteristic polynomial p(λ) are stable. Therefore, the equilibrium (0, 0) is a stable one. 5.2 Condition for the appearance of limit cycles Finally, it is noted that the appearance of the limit cycles in the FitzHugh–Nagumo neuron model is possible. For instance, if state feedback is implemented in the form I = −q x2 then (x1∗ , x2∗ ) = (0, 0) is again a fixed point for the model while the Jacobian matrix J at (x1∗ , x2∗ ) becomes   −γ  (27) J= −1 −3x22 + 2(1 + a)x2 − a − q The characteristic polynomial of the above Jacobian matrix is p(λ) = λ2 + (a + q + γ λ + (γ q + )). Setting the feedback gain q = −a − γ , the characteristic polynomial becomes p(λ) = λ2 +( 2 γ 2 +αγ −). Then, the condition (αγ − ) > 0 with a > 0,  > 0 and γ > 0 suffices for the appearance of imaginary eigenvalues. In such a case, the neuron model exhibits sustained oscillations (limit cycles).

6 Bifurcation analysis of the Neurospora circadian oscillator

The model’s parameters are defined as follows: vs denotes the rate of transcription of the frq gene into FRQ protein inside the nucleus. K i represents the threshold constant beyond which nuclear FRQ protein inhibits transcription of the frq gene. n is the Hill coefficient showing the degree of cooperativity of the inhibition process, vm is the maximum rate of frq mRNA degradation, and K M is the Michaelis constant related to the transcription process. Parameter K s defines the rate of FRQ synthesis and stands for the control input to the model. Parameter K 1 denotes the rate of transfer of FRQ protein from the cytoplasm to the nucleus, and K 2 denotes the rate of transfer of the FRQ protein from the nucleus into the cytoplasm. Parameter vd denotes the maximum rate of degradation of the FRQ protein, and K d is the Michaelis constant associated with this process. Fixed points are computed from the nullclines x˙1 = 0, x˙2 = 0, and x˙3 = 0. From Eq. (30) one has x3 =

K1 x2 K2

By substituting Eq. (31) and after intermediate operations, one obtains x1 also as a function of x2 that is x1 =

K 2 (vd + K 1 K d )x2 + K 1 K 2 K d x22 − (K 2 K d + K 2 x2 )K 1 x2 K 2 (K s K d + K s x2 )

By substituting Eqs. (31) and (32) into Eq. (28) and after intermediate operations, one obtains the equation −vm K 12 (kd − 1)x23 + [vm K in K 1 K 2 K d − vm K in K 1 K 2 +vm K 1 vd − vs K in K 1 K 2 K d + vs K in K 1 K 2 ]x22  + vm K in K 2 vd − vs K in K m K s K 2 − vs K 1n K 2 (vd + K 1 K d )  (33) + vs K in K 1 K 2 K d x2 − vs K in K m K 2 K s K d = 0 The control input of the circadian oscillator is taken to be zero, that is K s = 0. It is assumed that the bifurcating parameter is K d (Michaelis constant). Indicative values of the model’s parameters are selected according to Gonze et al. (2000). By substituting the numerical values of the parameters of the model into the characteristic polynomial, vs = 0.02 K i = 0.9, vm = 1.5, K m = 1.6, K s = 1, vd = 1.4, K 1 = 0.5, K 2 = 0.6, K c = 0.5 and setting n = 4, one obtains +(0.0039K d − 0.8196)} = 0

The model of the circadian cells is given by

123

(32)

x2 {0.3750(kd − 1)x22 + (0.2913K d + 0.7587)x2

6.1 Bifurcation diagrams of the fixed points

Kn x1 x˙1 = vs n i n − vm K i + x3 K m + x1 x2 x˙2 = −vd − K 1 x2 + K 2 x3 + x1 K s K d + x2 x˙3 = K 1 x2 − K 2 x3

(31)

(34)

A first fixed point is located at (28) (29) (30)

(x1∗,1 , x2∗,1 , x3∗,1 ) = (0, 0, 0)

(35)

For the binomial appearing in the above equation, the determinant is  = 0.0791K d2 + 1.6773K d − 0.6536 which

Biol Cybern (2014) 108:365–380

371

is positive for K d ≥0.3891. In such a case, there are two more fixed points at x2 =

−(0.2913K d + 0.7587) +

 0.0791K d2 + 1.6773K d − 0.6538

0.75(K d − 1)

(36)

x2 =

−(0.2913K d + 0.7587) −

 0.0791K d2 + 1.6773K d − 0.6538

0.75(K d − 1)

(37) The bifurcation diagram of the fixed points considering as bifurcation parameter the Michaelis constant K d is given in Fig. 1a. Next, the Jacobian matrix of the dynamic model of the system is computed, which is given by ⎞ ⎛ K in nx3n−1 KM 0 − (K n +x n )2 ⎟ ⎜−vm (K M +x1 )2 3 i ⎟ (38) vd K d J =⎜ 0 − − K1 K2 ⎠ ⎝ (K d +x2 )2

0

K1

−K 2

The associated characteristic polynomial is vd vm 2 p(λ) = λ3 + [ + K1 + K2 + ]λ Kd Km K 2 vd vm vd vm K 2 vd +[ + ( ) + K 1 + K 2 ]λ + . Kd Km Kd Km Kd (39) The variation ranges for parameter K d are taken to be K d ∈[K dmin , K dmax ] = [0.4, 0.9]. Then, the variation ranges of the rest of the parameters of the characteristic polynomial are as follows: c2 = [ Kvdd + K 1 + K 2 + Kvmm ]∈[3.59, 5.54], K 2 vd vm vd Kd + Km ( Kd ) K 2 vd K d ∈[0.88, 1.97].

c1 = vm Km

+ K 1 + K 2 ∈[3.42, 6.41] and c0 =

Since there is no change of sign in the first column of the Routh matrix, it can be concluded that the characteristic polynomial is a stable one. For polynomial p2 (λ), it holds λ3 | λ2 | λ1 | λ0 |

1 3.42 5.54 0.88 3.26 0 0.88

(43)

Since there is no change of sign in the first column of the Routh matrix, it can be concluded that the characteristic polynomial is a stable one. For polynomial p3 (λ), it holds λ3 | λ2 | λ1 | λ0 |

1 6.41 5.54 0.88 6.25 0 0.88

(44)

Since there is no change of sign in the first column of the Routh matrix, it can be concluded that the characteristic polynomial is a stable one. For polynomial p4 (λ), it holds λ3 | λ2 | λ1 | λ0 |

1 3.42 3.59 1.97 2.87 0 1.97

(45)

Since there is no change of sign in the first column of the Routh matrix, it can be concluded that the characteristic polynomial is a stable one. Consequently, the fixed point (x1∗,1 , x2∗,1 , x3∗,1 ) = (0, 0, 0) is a stable one (sink) for the considered variations of K d .

Next, the four Kharitonov polynomials are computed as follows:

6.2 Method to detect Hopf bifurcations in the circadian cells

p1 (λ) = c0max + c1max λ + c2min λ2 + c3min λ3

Using the characteristic polynomial of the linearized equivalent of the biological model, it is possible to formulate conditions for the emergence of Hopf bifurcations (Baird and Terman 2009). The whole procedure makes again use of the Routh–Hurwitz criterion. Consider the characteristic polynomial computed at a fixed point

p2 (λ) = c0min + c1min λ + c2max λ2 + c3max λ3 p3 (λ) = c0min + c1max λ + c2max λ2 + c3min λ3

(40)

p4 (λ) = c0max + c1min λ + c2min λ2 + c3max λ3 which take the values

P(λ) = cn λn + cn−1 λn−1 + cn−2 λn−2 + · · · + c1 λ + c0 (46)

p1 (λ) = 1.97 + 6.41λ + 3.59λ2 + λ3 p2 (λ) = 0.88 + 3.42λ + 5.54λ2 + λ3 p3 (λ) = 0.88 + 6.41λ + 5.54λ2 + λ3

(41)

p4 (λ) = 1.97 + 3.42λ + 3.59λ2 + λ3 The Routh criterion is applied for each characteristic polynomial pi (λ) i = 1, · · · , 4. For polynomial p1 (λ), it holds λ3 | 1 6.41 λ2 | 3.59 1.97 λ1 | 5.86 0 λ0 | 1.97

(42)

The coefficients ci , i = 0 · · · , n are functions of the bifurcating parameter, which is denoted as μ. The following matrix is introduced j = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝

cn−1 cn−3 cn−5 ···

cn

0

cn−2 cn−4 ···

cn−1 cn−3 ···

cn−(2 j−1) cn−(2 j−2) cn−(2 j−3)

⎞ ··· 0 ⎟ ··· 0 ⎟ ⎟ ··· 0 ⎟ ··· ··· ⎠ · · · cn−(2 j−n)

(47)

123

372

Biol Cybern (2014) 108:365–380

where coefficients cn−(2 j−m) with j = 1, 2, · · · , n and m = 1, 2, · · · , n become zero if n − (2 j − m) < 0 or if n − (2 j − m) > n. The Routh–Hurwitz determinants show the existence of eigenvalues which are conjugate imaginary numbers. A series of matrices containing the coefficients c j are defined as follows: 1 = cn−1   c c 2 = n−1 n cn−3 cn−2 ⎛ ⎞ cn−1 cn 0 3 = ⎝cn−3 cn−2 cn−1 ⎠ cn−5 cn−4 cn−3

vd vm 2 + K1 + K2 + ]λ Kd Km K 2 vd vm vd vm K 2 vd +[ + ( + K 1 + K 2 )]λ + Kd Km Kd Km Kd

p(λ) = λ3 + [

which after the substitution of numerical values for its parameters becomes     1.4 2.1525 0.7875 2 p(λ) = λ + +2.0375 λ + +1.0313 λ+ Kd Kd Kd (50) 3

Using the Routh–Hurwitz determinants, one has (48)

and so on. Each matrix  j , j = 1, . . . , n − 1 is square, and the first column contains every other coefficient of the characteristic polynomial cn−1 , cn−3 , · · · . The roots of the characteristic polynomial p(λ) have negative real parts if and only if det ( j ) > 0 for j = 1, . . . , n. The following two remarks can be made as follows: (i) Since det (n ) = c0 det (n−1 ), this means that a necessary condition for the roots of p(λ) to have negative real part is c0 > 0. If c0 = 0 then there is a zero eigenvalue. (ii) If det ( j ) > 0 for all j = 1, 2, · · · , n − 2 and if it holds det (n−1 ) = 0 then the characteristic polynomial contains imaginary roots. Provided that the rest of the roots of the characteristic polynomial have negative real part then Hopf bifurcations appear. An additional condition for the appearance of Hopf bifurcations is that the derivative of the determinant det (n−1 ) = 0 with respect to the bifurd det (n−1 ) = 0. cating parameter is not zero. That is dμ

1 = c1 =

According to the above criteria (i) and (ii), one can determine where possible saddle-node bifurcations (eigenvalue 0) and Hopf bifurcations (imaginary eigenvalues) appear. Therefore, one can formulate conditions about the values of the bifurcating parameter that result in a particular type of bifurcations. Equivalently, conditions about the appearance of Hopf bifurcations can be obtained from the Routh table according to Example 4, presented in Sect. 3.1. The zeroing of certain rows of the Routh table gives an indication that the characteristic polynomial contains imaginary conjugate eigenvalues. By requiring the elements of these rows of the Routh matrix to be set equal to zero, one finds the values of the bifurcating parameter that result in the appearance of Hopf bifurcations. For the previous characteristic polynomial of the circadian cells, one has

1.4 + 2.0375 Kd

(51)

For 1 > 0, it should hold that K d < −0.6871 which is not possible since K d is a positive parameter. For the determinant of 2 to be zero, where    1.4  1 c1 c0 K d + 2.0375 2 = = (52) 0.7875 2.1525 c3 c2 Kd K d + 1.0313 one obtains the relation

3.0135 K d2

+

5.0420 Kd

+ 2.1013 = 0 which

in turn gives K d = −1.2722, K d = −1.1273. Again, these values are not acceptable since K d > 0. Therefore, no Hopf bifurcations can appear at the fixed point (x1∗,1 , x2∗,1 , x3∗,1 ) = (0, 0, 0), under the existing set of parameters’ values. At the same conclusion, one arrives if the Routh–Hurwitz table is used, that is λ3 | λ2 | λ1 | λ0 |

123

(49)

1 3.0135 K d2

1.4 Kd

−

+ 2.0375 5.0420 K d − 2.1013 0.7875 Kd

2.1525 Kd 0.7875 Kd

0

(53)

For the first coefficient of the second row of the matrix, it holds 1.4 K d +2.0375 > 0. To obtain zero coefficients at the line associated with λ1 (so as imaginary eigenvalues to appear in the characteristic polynomial of the preceding line), one has + 5.0420 again the relation 3.0135 K d + 2.1013 = 0 which in turn K2 d

gives K d = −1.2722, K d = −1.1273. As mentioned, these values are not acceptable since K d > 0. Therefore, no Hopf bifurcations can appear at the fixed point (x1∗,1 , x2∗,1 , x3∗,1 ) = (0, 0, 0), under the existing set of parameters’ values.

7 Feedback control of bifurcations Next, it is considered that input K s is a feedback control term of the form K s = c2 x2 . Then, the dynamics of the circadian

Biol Cybern (2014) 108:365–380

373

oscillator becomes x˙1 = x˙2 =

3rd fixed point:

Kn x1 vs n i n − vm K i + x3 K m + x1 x2 −vd − K 1 x2 + K 2 x3 K d + x2

+ c2 x1 x3

(54)

x˙3 = K 1 x2 − K 2 x3

Feedback can change the oscillation characteristics of the circadian model. It can change the type of stability characterizing the system’s fixed points. Additionally, fixed-point bifurcations can be altered through state feedback control. The previously analyzed Eq. (33) is considered again for computing fixed points  −vm k12 (K d − 1)x23 + vm K in K 1 K 2 K d − vm K in K 1 K 2  +vm K 1 vd − vs K in K 1 K 2 K d + vs K in K 1 K 2 x22  + vm K in K 2 vd − vs K in K m K s K 2 − vs K 1n K 2 (vd + K 1 K d )  (55) + vs K in K 1 K 2 K d x2 − vs K in K m K 2 K s K d = 0 Next, the nominal value of K d will be considered known and will be substituted in the characteristic polynomial, while the input K s will be considered equal to the feedback law of the form K s = c2 x2 . In such a case, the equation for computing x2 becomes −vm k12 (K d − 1)x23 + [vm K in K 1 K 2 K d − vm K in K 1 K 2 +vm K 1 vd − vs K in K 1 K 2 K d + vs K in K 1 K 2 −vs K in K m c2 K 2 ]x22 + [vm K in K 2 vd −vs K 1n K 2 (vd + K 1 K d ) + vs K in K 1 K 2 K d −vs K in K m K 2 c2 K d ]x2 = 0

(56)

Substituting numerical values for the model’s parameters and taking K d = 1.4, the equation for computing x2 becomes x2 {−0.15x22 + (1.1665 + 0.0126c2 )x2 +(0.8157 + 0.0176c2 )} = 0

(57)

The fixed points for the model of the circadian cells are found to be: 1st fixed point: (x1∗,1 , x2∗,1 , x3∗,1 ) = (0, 0, 0), whereas from the computation of the associated determinant, one obtains  = ((1.6665+0.0126c2 )2 +0.6(0.8157+0.0176c2 )) which is positive if c2 takes values in the intervals c2 < −251.15 and c2 > −82.39. Thus, for  > 0, one has the fixed points expressed as functions of the feedback gain c2 . 2nd fixed point: x1∗,2 =

K 1 (K d − 1)x2 ∗,2 + vd c2 (kd + 1)

x2∗,2 =

−(1.665 + 0.0126c2 ) +

x3∗,2 =



1.58·10−4 c22 + 0.0527c2 + 3.2694 −0.3

K 1 ∗,2 x2 K2

K 1 (K d − 1)x2 ∗,3 + vd c2 (kd + 1)

x1∗,3 =

(58)

−(1.665 + 0.0126c2 ) −

x2∗,3 =

 1.58·10−4 c22 + 0.0527c2 + 3.2694 −0.3

K 1 ∗,3 x2 x3∗,3 = K2

(59)

The bifurcation diagram of the fixed points considering as bifurcating parameter the feedback control gain c2 is given in Fig. 1b. Next, from the right part of Eq. (54), the system’s Jacobian is computed J=

⎛ Km ⎜−vm (K +x )2 m 1 ⎜ ⎜ c2 x2 −vd ⎝

0

0 Kd − K 1 + c2 x1 (K d +x2 )2 K1

2⎞ −K n nx n−1 vs (Kin +x3n ) ⎟ ⎟ 3 i

K2

⎟ ⎠

(60)

−K 2

From the system’s Jacobian, a generic form of the characteristic polynomial can be computed at the fixed point (x1∗ , x2∗ , x2∗ ). After intermediate operations, one obtains det (λI − J ) = vd kd Km K 1 − c2 x1 ] + [vm + ]}λ2 λ3 + {[k2 + 2 K d + x2 K m + x1 2  vd kd Km + [ K − K 1 K 2 ] + [vm + ]· 2 1 (K d + x2 ) (K m + x1 )2  vs K in nx3n−1 vd K d [K 2 + ( K 1 − c2 x1 )] − λ (K in + x3n )2 (K d + x2 )2  Km vd K d ][( K − c2 x1 ) − K 1 K 2 ] + [vm + 2 2 1 (K K m + x1 d + x2 )  vs K in nx3n−1 vd K d − ( ) + K 1 − c2 x2 (61) (K in + x3n )2 K d + x2 2 or equivalently det (λI − J ) = λ3 + c f2 λ2 + c f1 λ + c f0

(62)

where   vd K d − K K − c x K 1 2 1 1 2 (K d + x2 )2 K m + x1 2   n−1  n vs K i nx3 vd K d − n + K 1 − c2 x2 (K i + x3n )2 K d + x2 2     vd kd Km + v · c f1 = K − K K + m 1 1 2 (K d + x2 )2 (K m + x1 )2     vs K in nx3n−1 vd K d − · K2 + K − c x 1 2 1 (K in + x3n )2 (K d + x2 )2     vd kd Km + v c f2 = k2 + K − c x + m 1 2 1 K d + x2 2 K m + x1 2 

c f0 =

vm +

Km

 

(63)

(64) (65)

123

374

Biol Cybern (2014) 108:365–380 25

The associated Routh matrix is 1 c f1 λ3 | λ2 | c f2 c f0 λ1 | c f 0 − c f 1 c f 2 0 c f0 λ0 |

(66)

p1 (λ) = 23.67 + 70.39λ + 31.31λ + λ 2

p2 (λ) = 23.27 + 68.33λ + 31.91λ + λ 2

15

x2

The condition for obtaining Hopf bifurcation is c f2 > 0 and c f0 − c f1 c f2 = 0. It can be confirmed that for the fixed point (x1∗,2 , x2∗,2 , x3∗,2 ) and for the previously used set of values of the parameters of the circadian oscillator model, solutions c2 obtained from the condition c f0 − c f1 c f2 = 0 result in c f2 < 0; therefore, Hopf bifurcations do not appear. The bifurcation diagram of Fig. 1b shows the limit values of c2 beyond which x2∗,2 remains positive. Such a value is c2 < −45.2. It is examined next; if for specific range of variations of parameter c2 , the fixed point (x1∗,2 , x2∗,2 , x3∗,2 ) is a stable or an unstable one. To this end, Kharitonov’s theorem will be used again. For variation range c2 ∈[−60.2, −50.2], it holds that c f2 ∈[31.31, 31.91], c f1 ∈[68.33, 70.39] and c f0 ∈[23.27, 23.67]. The associated four Kharitonov polynomials take the values

20

10

5

0

0.4

0.5

0.6

0.7

0.8

0.9

1

Kd

(a)

3 14

3

p3 (λ) = 23.27 + 70.39λ + 31.91λ2 + λ3

(67)

12

p4 (λ) = 23.67 + 68.33λ + 31.31λ2 + λ3

10

The Routh table for Kharitonov polynomial p1 (λ) is

(68)

2

8

x

λ3 | 1 70.39 λ2 | 31.31 23.67 λ1 | 69.63 0 λ0 | 23.67

6

Since there is no change of sign of the coefficients of the first column of the Routh table, it can be concluded that the characteristic polynomial p1 (λ) is a stable one. The associated Routh table for the second Kharitonov polynomial p2 (λ) is

4

λ3 | 1 68.33 λ2 | 31.91 23.27 λ1 | 67.59 0 λ0 | 23.27

0

(69)

−80

−60

−40

−20

0

20

(b) Fig. 1 a Bifurcation of fixed point x2∗,2 (blue line) and x2∗,3 (red line) of the circadian oscillators under variation of the Michaelis parameter K d , b Bifurcation of fixed point x2∗,2 (blue line) and x2∗,3 (red line) of the circadian oscillators under state feedback with gain c2

(70)

Since there is no change of sign of the coefficients of the first column of the Routh table, it can be concluded that the

123

−2 −100

c2

Since there are no changes of sign in the coefficients of the first column of the Routh tables, it can be concluded that polynomial p2 (λ) is also stable. The associated Routh table for Kharitonov polynomials p3 (λ) and p4 (λ) are λ3 | 1 70.39 λ2 | 31.91 23.27 λ1 | 69.66 0 λ0 | 23.27

2

characteristic polynomial p3 (λ) is a stable one. Finally, the associated Routh table for the fourth Kharitonov polynomial

Biol Cybern (2014) 108:365–380

375

p4 (λ) is

2

λ3 | 1 68.33 λ2 | 31.31 23.67 λ1 | 67.57 0 λ0 | 23.67

1

(71)

x −x

* 1

1.5

1

0.5

0

0

5

10

p2 (λ) = −88.17 − 175.03λ − 55.12λ + λ

20

* 2 2

1

0

−1

0

5

10

time

(a) 3

(72)

p4 (λ) = −84.13 − 175.03λ − 58.12λ2 + λ3

2.5

The associated Routh table for the first Kharitonov polynomial p1 (λ) is * 2

2

1.5

x 2 −x

(73)

1

Since there is change of sign of the coefficients of the first column of the Routh table, it can be concluded that the characteristic polynomial p1 (λ) is an unstable one. The associated Routh table for Kharitonov polynomial p1 (λ) and p2 (λ) are | 1 −175.03 λ2 | −55.12 −88.17 0 λ1 | −176.63 λ0 | −88.17

15

2

3

p3 (λ) = −88.17 − 163.82λ − 55.12λ2 + λ3

λ3 | 1 −163.82 λ2 | −58.52 −84.14 0 λ1 | −165.25 λ0 | −84.14

20

3

p1 (λ) = −84.13 − 163.82λ − 58.52λ2 + λ3 2

15

time

x −x

Since there are no changes of sign in the coefficients of the first column of the Routh tables, it can be concluded that polynomial p4 (λ) is stable. From the previous analysis, it can be seen that all Kharitonov polynomials have stable roots; therefore, the fixed point (x1∗,2 , x2∗,2 , x3∗,2 ) of the circadian oscillator under feedback gain c2 ∈[−60.2, −50.2] is a stable one. Next, the case that the feedback gain c2 ∈ [−270.2,−260.2] is examined. The variation ranges for the coefficients of the characteristic polynomial are c f2 ∈ [−58.52, −55.12], c f1 ∈[−175.03, −163.82] and c f0 ∈ [−88.17, −84.13]. The four Kharitonov polynomials of the system are as follows:

0.5

0

λ3

(74)

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

*

x −x 1

From the changes of sign of the coefficients of the first column of the Routh tables, it can be concluded that polynomial p2 (λ) is unstable. The associated Routh table for the third Kharitonov polynomial p3 (λ) is λ3 | 1 −163.82 2 λ | −55.12 −88.17 0 λ1 | −165.42 λ0 | −88.17

−0.5

(75)

Since there is change of sign of the coefficients of the first column of the Routh table, it can be concluded that the characteristic polynomial p3 (λ) is an unstable one. Finally, the

1

(b) Fig. 2 Fixed point of the FitzHugh–Nagumo neuron is a stable one: a diagrams of the state variables, b phase diagrams

associated Routh table for the fourth Kharitonov polynomial p4 (λ) is 1 −175.03 λ3 | λ2 | −58.12 −84.13 0 λ1 | −176.48 λ0 | −84.13

(76)

123

376

x −x * 1 1

20

10

0

0

5

10

15

20

15

20

15

20

time 20

8 Simulation tests

2

10

0

0

5

10

time

3

x −x *

3

20

10

0

0

5

10

time

(a) 20 15

x −x * 2 2

Remark 1 Either by applying the procedure that is explained in Sect. 6.2 and which is based on the computation of the Routh–Hurwitz determinants or by finding conditions for zeroing rows of the Routh table as explained again in Sect. 6.2, one can compute the values of the uncertain parameters of the system which result in Hopf bifurcations at the fixed point (which means that the dynamic behavior of the system is characterized by sustained oscillations). As explained in the example given in Sect. 6.2, the Hopf bifurcations appear for specific values of the uncertain parameters. In the relevant bibliography, one can find results on the appearance of such an oscillatory behavior in the dynamics of the circadian oscillators (François 2005; Gonze et al. 2002; Leloup et al. 1999; Zhang et al. 2010). However, when moving out of these values, the oscillations cease to exist.

x −x *

2

From the changes of sign of the coefficients of the first column of the Routh tables, it can be concluded that the characteristic polynomial p4 (λ) is unstable. From the previous analysis, it can be seen that all Kharitonov polynomials have unstable roots; therefore, the fixed point (x1∗,2 , x2∗,2 , x3∗,2 ) of the circadian oscillator under feedback gain c2 ∈[−270.2, −260.2] is an unstable one.

Biol Cybern (2014) 108:365–380

10 5

123

0

0

2

4

6

8

10

12

*

x −x 1

1

15

x −x * 3 3

The following example is concerned with bifurcations of the fixed point of the FitzHugh–Nagumo neuron. Assuming zero external current input and known parameters of the model, it has been found that the associated fixed point is (x1∗,1 , x2∗,1 , x2∗,1 ) = (0, 0, 0). As confirmed by the results given in Fig. 2, this fixed point is a stable one. The nominal values of the model’s parameters were  = 1.5, γ = 0.1 and a = 7. Next, examples on the stability of fixed points in the circadian oscillator model are given under variation of specific parameters of the dynamical model. The nominal values of the model’s parameters were as follows: vs = 0.02, K i = 0.9, vm = 1.5, K m = 1.6, vd = 1.4 K 1 = 0.5, K 2 = 0.6, n = 4, while K d was an uncertain parameter. The control input of the circadian oscillator is taken to be zero, that is K s = 0. It was assumed that the Michaelis constant parameter varies within intervals, which was K d ∈[0.4, 0.9], and based on this interval, the ranges of variation of the Jacobian matrix of the model computed at the fixed point (x1∗,1 , x2∗,1 , x2∗,1 ) = (0, 0, 0) were found. By applying Kharitonov’s theorem, it was possible to show that conditions about the stability of the fixed point were satisfied. The associated results are depicted in Fig. 3. Moreover, the case of state feedback was considered in the form K s = c2 x2 and using parameter c2 as feedback gain, while it was taken that K d = 1.4. The variation

10

5

0

0

5

10

15

20

*

x −x 2

2

(b) Fig. 3 First fixed point of the circadian oscillator becomes a stable one for parameter K d ∈[0.4, 0.9]: a diagrams of the state variables, b phase diagrams

range of c2 was taken to be the interval [−60.2, −50.2], and based on this, the ranges of variation of the coefficient of the characteristic polynomial of the system’s Jacobian matrix were computed. Stability analysis followed, using Kharitonov’s theorem. First, it was shown that by choosing

Biol Cybern (2014) 108:365–380

377

x 1 −x *1

1.4

*

x 1 −x 1

2

1

0

0

5

10

15

1.2

1

20

0

0.02

0.04

0.06

0.12

0.14

0.16

0.1

0.12

0.14

0.16

0.1

0.12

0.14

0.16

1.25

1.3

1.35

1.4

2.5

2.6

2.7

2.8

2

x 2 −x *

*

x 2 −x 2

0.1

3

5

0

−5

2.5

2

0

5

10

15

20

0

0.02

0.04

0.06

0.08

time

time 6

*

* x −x

3

5

x −x 3 3

0.08

time

time

−5

3

0

4

2

0

5

10

15

20

0

0.02

0.04

0.06

0.08

time

time

(a)

(a) 2.8

3

2.6

2

1

2.4

2

x −x *

2

x −x *

2

2

2.2

0 −1

2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

1

1.05

1.1

1.15

x −x

x −x 1

1.2 *

*

1

1

1

6

2

* 3

0

x −x

5

3

*

x −x 3 3

4

4

−2 3

−4 −0.5

0

0.5

1

1.5

2

2.5

*

x −x 2

2

(b) Fig. 4 Second fixed point of the circadian oscillator is stable for feedback gain value c2 ∈[−60.2, −50.2]: a diagrams of the state variables, b phase diagrams

−60.2 < c2 < −50.2, e.g., c2 = −55; then, the fixed point (x1∗,2 , x2∗,2 , x3∗,2 ) remained stable. Next, it was shown that by setting −270.2 < c2 < −260.2, e.g., c2 = −265, the fixed point became unstable. Indeed in the first case, that is when c2 = −55, the roots of the characteristic polynomial that is associated with the

2

2.1

2.2

2.3

2.4 *

x −x 2

2

(b) Fig. 5 Second fixed point of the circadian oscillator is unstable one for feedback gain value c2 ∈[−270.2, −260.2]: a diagrams of the state variables, b phase diagrams

Jacobian of the system are λ1 = −29.49, λ2 = −1.52 and λ3 = −0.42, which means that the fixed point is a stable one. In the second case, that is for c2 = −265, the roots of the characteristics polynomial became λ1 = 59.5, λ2 = −2.21 and λ3 = −0.65 which means that the fixed point is an unstable one. The associated results are depicted in Figs. 4 and 5.

123

378

Biol Cybern (2014) 108:365–380

9 Conclusions

10.1 The Routh–Hurwitz stability criterion

The paper has studied bifurcations in biological models, such as circadian cells. It has been shown that the stages for analyzing the stability features of the fixed points that lie on bifurcation branches comprise (i) the computation of fixed points as functions of the bifurcation parameter and (ii) the evaluation of the type of stability for each fixed point through the computation of the eigenvalues of the Jacobian matrix that is associated with the system’s nonlinear dynamics model. Stage (ii) requires the computation of the roots of the characteristic polynomial of the Jacobian matrix. This problem is nontrivial since the coefficients of the characteristic polynomial are functions of the bifurcation parameter and the latter varies within intervals. To obtain a clear view about the values of the roots of the characteristic polynomial and about the stability features they provide to the system, the use of interval polynomials theory and particularly of Kharitonov’s stability theorem has been proposed. In this approach, the study of the stability of a characteristic polynomial with coefficients that vary in intervals is equivalent to the study of the stability of four polynomials with crisp coefficients computed from the boundaries of the aforementioned intervals. The method for fixed-point stability analysis was tested on two different biological models: (i) the FitzHugh–Nagumo neuron model, (ii) a model of circadian cells which produce proteins out of an RNA transcription and translation procedure. The efficiency of the proposed stability analysis method has been confirmed through numerical and simulation tests.

To evaluate the stability of the roots of the characteristic polynomial

Appendix 10 A Stability analysis through the Routh–Hurwitz criterion Since the Routh–Hurwitz criterion is meanwhile of the category well known but not widely known, a quick look at the ideas will be taken behind this classical method and a few concrete examples will be analyzed as illustrations of how to apply it. The Routh–Hurwitz criterion provides a systematic method for finding out if all the roots of a polynomial lie in the left complex semi-plane. In the case of a dynamical system described in the s-frequency domain by a transfer function, one can find if the denominator’s characteristic polynomial contains only stable roots thus concluding that the system is also stable. Moreover, one can find out the number of complex eigenvalues which lie in the right complex semi-plane and thus are unstable. Additionally, one can find if there are roots of the polynomial with purely imaginary part, resulting a sustained oscillations dynamics (limit cycles).

123

p(λ) = an λn + an−1 λn−1 + · · · + a1 λ + a0

(77)

the following Routh table is considered: s n | an s n−1 | an−1 s n−2 | q1 (n−2) s n−3 | q1 (n−3) ··· | ··· s 1 | q1(1) (0) s 0 | q1

an−2 an−3

q2 (n−2) q2 (n−3) ··· q2(1)

an−4 · · · an−5 · · · q3 (n−2) · · · q3 (n−3) · · · ··· ···

(78)

where  (n−2)

q1

(n−3)

q1

det

=−

=−

⎛

an an−2 an−1 an−3 an−1



 (n−2)

⎞

an−1 an−3 det ⎝ (n−2) (n−2) ⎠ q1 q2 q1 (n−2)

q2

(n−3)

q2

=

det

=−

an an−4 an−1 an−4



an−1 ⎛ ⎞ an−1 an−4 (79) det ⎝ (n−2) (n−2) ⎠ q1 q3 − q1 (n−2)

etc. The Routh–Hurwitz stability analysis points out that the number of unstable roots of the characteristic polynomial (roots with positive real part) is equal to the number of changes of sign appearing in the first column of the Routh table. Example 1 The characteristic polynomial p(s) = s 4 +5s 3 + 20s 2 + 30s + 40 is considered. The associated Routh table is s4 | 1 s3 | 5 s 2 | 14 s 1 | 15.714 s 0 | 40

20 40 30 0 40 0 0

(80)

No change of sign is observed in the first column of the Routh table. Therefore, the characteristic polynomial has no roots in the right complex semi-plane and is a stable one. Example 2 The characteristic polynomial p(s) = s 3 +2s 2 + 4s + 30 is considered. The associated Routh table is s3 | 1 4 s 2 | 2 30 s 1 | −11 0 s 0 | 30

(81)

There are two changes if sign in the first column of the Routh table; therefore, the characteristic polynomial has two unstable roots.

Biol Cybern (2014) 108:365–380

379

Example 3 The characteristic polynomial p(s) = s 5 +2s 4 + 3s 3 + 6s 2 + 12s + 18 is considered. In this case, an element in the first column of the Routh table becomes equal to 0. This is substituted by a small positive number . Next, the sign of the elements of the first column of the Routh table is examined while →0. Thus, the associated Routh table becomes s5 s4 s3 s2 s1 s0

| | | | | |

1 2  6−

6  18−18−18 2 6−6

3 6 3 18 0

12 18 0 0 0

(82)

18

For →0+ , the elements of the first column become 1, 2,  > 0, 6 − 6 < 0, 3, 18. There are two changes of sign in the elements of the first column of Routh table, and consequently, there are two unstable roots in the right complex semi-plane. Example 4 The characteristic polynomial p(s) = s 5 +6s 4 + 10s 3 + 35s 2 + 24s + 44 is considered. In such a case, complete row of the associated Routh table can become equal to zero. This happens when among the roots of the characteristic polynomial there are conjugate imaginary ones or complex roots of opposite sign. To process this case, an auxiliary characteristic polynomial is generated using the row that precedes the row where zeroing appeared. By dividing the initial characteristic polynomial with the auxiliary characteristic polynomial, one obtains a new characteristic polynomial for which the computation of the Routh table is repeated. The Routh table of the initial characteristic polynomial is s5 s4 s3 s2 s1 s0

| | | | | |

1 10 24 6 35 44 25 50 6 3 0 11 44 0 0 0 0

(83)

The auxiliary characteristic polynomial is p(s) = 11s 2 + 44 with purely imaginary roots s1,2 = ± j2. The new characteristic polynomial is computed using q(s) = p(s)/a(s) = s 3 + 6s 2 + 6s + 11. For q(s), the Routh table is s3 s2 s1 s0

| 1 6 | 6 11 | 25 6 0 | 11

(84)

For the new Routh table, there is no change of sign in its first column; thus, all roots of q(s) lie in the left complex semi-plane. Consequently, for the initial characteristic

polynomial, it holds p(s) = (11s 2 + 44)q(s) which means that p(s) has two purely imaginary roots (Hopf bifurcations) while the rest of its roots are stable. Example 5 Consider the characteristic polynomial p(s) = s 6 + 2s 5 + 9s 2 + 12s 3 + 43s 2 + 50s + 75. The associated Routh table is s6 | 1 s5 | 2 s4 | 3 s3 | 0 s2 | s1 | s0 |

9 12 18 0

43 50 75 0

75 0 0 0

(85)

There is zeroing of one row. The auxiliary characteristic polynomial is formulated using the row preceding to zeroing, that is a(s) = 3s 4 + 18s 3 + 75 = 3(s 4 + 6s 2 + 25). By diving the initial characteristic polynomial with the auxiliary characteristic polynomial, the new characteristic polynomial is obtained. That is 1 p(s) = (s 2 + 2s + 3)⇒ a(s) 3 1 2 p(s) = (s + 2s + 3)(s 4 + 6s 2 + 25) 3

q(s) =

(86)

The polynomial q1 (s) = s 4 + 6s 2 + 25 has complex roots (every two of them are opposite to each other), that is ρ1 = 1 + j2, ρ2 = 1 − j2, ρ3 = −1 + j2 and ρ4 = −1 − j2. The = s 2 + 2s + 3 has the polynomial q2 (s) √ √ complex conjugate roots ρ1 = −1+ j 2 and ρ2 = −1− j 2. From the previous analysis, it can be seen that the characteristic polynomial p(s) has roots in the right complex semi-plane and is unstable. 10.2 Use of the Routh–Hurwitz criterion for the computation of relative stability By performing the change of variable s→ p − a in the characteristic polynomial, one can find the relative position of the roots of the characteristic polynomial with respect to a reference value. This means that one can conclude if the roots of the characteristic polynomial are located to the left of the reference eigenvalue a (they decay faster than a) or to the right of the reference eigenvalue (they decay slower than a). The following characteristic polynomial is considered: p(s) = s 3 + 10.2s 2 + 21s + 2 and the associated Routh table becomes s 3 | 1 21 s 2 | 10.2 2 s 1 | 20.804 0 s0 | 2

(87)

123

380

Biol Cybern (2014) 108:365–380

The system has the following stable poles p1 = −7.39, p2 = −2.70, p3 = −0.10. It is examined if the poles of the system decay faster than the eigenvalue, a = −0.2. The change of variables s = p − 0.2 or p = s + 0.2 is considered. The new characteristic polynomial is q(s) = p 3 + 9.6 p 2 + 16.32 p − 1.728 and the associated Routh table becomes p3 p2 p1 p0

| 1 16.32 | 9.6 −1.728 | 16.580 0 | −1.732

(88)

It can be noticed that there is one change of sign in the first column of the Routh table; therefore, the characteristic polynomial has an eigenvalue which lies to the right of the reference eigenvalue a = −0.2.

References Aqil M, Hong KS, Jeong MY (2012) Synchronization of coupled chaotic FitzHugh–Nagumo systems. Commun Nonlinear Sci Numer Simul 17:1615–1627 Baird Emertrout G, Terman DH (2009) Series in interdisciplinary applied Mathematics. Mathematical foundations of neuroscience. Springer, Berlin Chen G, Moiola JL, Wang HO (2000) Bifurcation control: theories methods and applications. Int J Bifurc Chaos 10(3):511–548 Clairambault J (2006) Physiologically based modelling of circadian control on cell proliferation. In: 28th IEEE EMBS annual international conference, New York, USA, Sep 2006 Clairambault J (2008) A step toward optimization of cancer therapeutics: physiologically based modeling of circadian control on cell proliferation. IEEE Eng Med Biol Mag 27(1):20–24 Dorak RO (2013) Control of repetitive firing in Hodgkin–Huxley nerve fibers using electric fields. Chaos, Solitons and Fractals 52:66–72 Duan H, Cai C, Han C, Che Y (2012) Bifurcation control in Morris– Lecar neuron model via washout filter with a linear term based on filter-aided dynamic feedback. Adv Mater Res 485:600–603 Efimov DV, Fradkov AL (2006) Adaptive tuning to bifurcation for timevarying nonlinear systems. Automatica 42:417–425 Francois P (2005) A model for the Neurospora circadian clock. Biophys J 88:2369–2383 Gerstner W, Kistler W (2002) Spiking neuron models: single neurons, populations, plasticity. Cambridge University Press, Cambridge Gonze D, Leloup JC, Goldbeter A (2000) Theoretical models for circadian rhythms in Neurospora and Drosophila. C R Acad Sci Paris, Sciences de la Vie / Life Sciences 323:57–67 Gonze D, Halloy José (2002) Robustness of circadian rhythms with respect to molecular noise. Proc Natl Acad Sci 99(2):673–678 Gonze D, Halloy J, Leloup JC, Goldbeter A (2003) Stochastic models for circadian rhythmes: effect of molecular noise on periodic and chaotic behaviour. C R Biol 326:189–203 Haschke R, Steil JJ (2005) Input-space bifurcation manifolds of recurrent neural networks. Neurocomputing 64:25–38 Lysetskiy M, Zurada JM (2004) Bifurcating neuron: computation and learning. Neural Netw 17:225–232 Nagy B (2008) Comparison of the bifurcation curves of a two-variable and a three-variable circadian rhythm model. Appl Math Model 38:1587–1598

123

Laroche B, Claude D (2004) Flatness-based control of PER protein oscillations in a Drosophila model. IEEE Trans Autom Control 49(2):175–183 Leloup JC, Gonze D, Goldbeter A (1999) Limit cycle models for circadian rhythms based on transcriptional regulation in Drosophila and Neurospora. J Biol Rhythms 14(6):433–448 Leloup JC, Gonze D, Goldbeter A (2006) Computational models for circadian rythms: deterministic versus stochastic approaches. Computational Systems Biology. Elsevier, Amsterdam Lee JH, Sancar A (2011) Circadian clock disruption improves the efficacy of chemotherapy through p73-mediated apoptosis. In: Proceedings of the National Academy of Sciences Lévi FA (2008) The circadian timing system: a coordinator of life processes. IEEE Eng Med Biol Mag 27(1):17–19 Liao X, Wong KW, Wu Z (2001) Bifurcation analysis of a two-neuron system with distributed delays. Phys D 149:123–141 Nagy B (2009) Analysis of the biological clock of Neurospora. J Comput Appl Math 226:298–305 Nguyen LH, Hong KS (2012) Hopf bifurcation control via a dynamic state-feedback control. Phys Lett A 376:442–446 Nguyen LH, Hong KS, Park S (2012) Bifurcation control of the Morris– Lecar neuron model via a dynamic state feedback ontrol. Biol Cybern 106:587–594 Rigatos GG, Siano P (2011) Design of robust electric power system stabilizers using Kharitonov’s theorem. Math Comput Simul 82(1):181–191 Rigatos G, Rigatou E (2013) A Kalman Filtering approach to robust synchronization of coupled neural oscillators. In: ICNAAM 2013, 11th international conference on numerical analysis and applied mathematics, Rhodes, Greece, Sep. 2013 Rigatos GG (2013) Advanced models of neural networks: nonlinear dynamics and stochasticity in biological neurons. Springer, Berlin Song Y, Han M, Wei J (2005) Stability and Hopf bifurcation analysis on a simplified BAM network with delays. Phys D 200:185–204 That LT, Ding Z (2012) Reduced-order observer design of multi-output nonlinear systems with application to a circadian model. Trans Inst Meas Control 35(4):417–425 Toscano R, Lyonnet P (2010) Robust static output feedback controller synthesis using Kharitonov’s theorem and evolutionary algorithms. Inf Sci 180:20232028 Tsumoto K, Yoshinaga T, Iida H, Kawakami H, Ashara K (2006) Bifurcations in a mathematical model for circadian oscillations of clock genes. J Theor Biol 239:101–122 Vidal A, Zhang Q, Médigue C, Fabre S, Clément F (2012) DynPeak: An Algorithm for Pulse Detection and Frequency Analysis in Hormonal Time Series. PLoS One 7(7):e39001 Wang J, Chen L, Fei X (2004) Bifurcation control of the HodgkinHuxley equations. Chaos, Solitons and Fractals 33:217–224 Wang H, Yu Y, Zhu R, Wang S (2013) Two-parameter bifurcation in a two-dimensional simplified Hodgkin–Huxley model. Commun Nonlinear Sci Numer Simul 18:184–193 Xu X, Hu HY, Wang HL (2006) Stability switches, Hopf bifurcations and chaos of a neuron model with delay dependent parameters. Phys Lett A 354:126–136 Zhang J, Bierman A, Wen JT, Julius A, Figueiro M (2010) Circadian system modeling and phase control. In: 49th IEEE conference on decision and control, Dec. 2010 Zhang T, Wang J, Fei X, Deng B (2007) Synchronization of coupled FitzHugh–Nagumo systems via MIMO feedback linearization control. Chaos Solitons and Fractals 33(1):194–202