Stability analysis for continuous-time and discrete-time genetic regulatory networks with delays

Applied Mathematics and Computation 274 (2016) 628–643

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Stability analysis for continuous-time and discrete-time genetic regulatory networks with delays Tingting Liu, Xian Zhang, Xiangyu Gao∗ School of Mathematical Science, Heilongjiang University, Harbin 150080, China

a r t i c l e

i n f o

Keywords: GRNs Exponential stability Dini derivative method Semi-discretization technique IMEX-θ method Linear matrix inequality(LMI) method

a b s t r a c t This paper investigates the global exponential stability problem for the genetic regulatory networks (GRNs) with delays in both the continuous-time case and the discrete-time case. First, Dini derivative method is applied for a new Lyapunov functional to obtain a sufficient condition of the global exponential stability for a class of the continuous-time GRNs, which is given in the form of several elementary inequalities. This sufficient condition is simple to be easily implemented on the computer. Second, by using the semi-discretization technique and the IMEX-θ method, two discrete-time GRN systems have been derived, which can preserve the dynamical characteristics of the continuous-time systems. Furthermore, it is shown that under the same sufficient conditions obtained earlier, these two discrete-time GRN systems are globally exponentially stable. Finally, a pair of examples are given to show the validity of the obtained stability analysis results. © 2015 Elsevier Inc. All rights reserved.

1. Introduction During the past two decades, GRNs have become an important new area of research in the biological and biomedical sciences, a large amount of outstanding results have been published in recent years (see [1–4] and the references therein). Different kinds of mathematical models have been applied to investigate the behaviors of GRNs, for example, Bayesian network models [5], Petri net models [6], Boolean models [7], the discrete model and differential equation models [8]. Among these models, the differential equation model describes the changing rate of the concentrations of gene products, such as mRNAs and proteins, as continuous values. The discrete-time systems are very important in the practical experiment implementations or in the computational simulations of the continuous-time GRNs. Time delay should be considered in GRNs due primarily to the slow processes of transcription, translation, translocation as well as the finite switching speed of amplifiers. It has been pointed out in [8–11] that delays may play an important role in dynamics of genetic networks, and mathematical models without addressing the delay effects may even provide wrong predictions of the mRNAs and proteins concentrations [8,10,12]. It has also been shown in [13] that, by mathematically modeling recent data, the observed oscillatory expression and activity of three proteins are most likely driven by transcriptional delays. On the other hand, stability is one of the most essential properties for designing or controlling GRNs [1], since stability is related to the ability of an organism to robustly regulate its function in spite of moving the state of the organism away from equilibrium points. An important issue in GRNs is to analyze stability of equilibrium points for the continuous-time and discretetime models. Due to this complexity of GRNs, the problem of stability analysis for GRNs has recently stirred increasing research ∗

Corresponding author. Tel.: +86 45186608282; fax: +86 45186608240. E-mail address: [email protected], [email protected], [email protected] (X. Gao).

http://dx.doi.org/10.1016/j.amc.2015.11.040 0096-3003/© 2015 Elsevier Inc. All rights reserved.

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

629

interest [9,14–19]. The stability criteria of GRNs with time delay are divided into two types: the delay-dependent stability criteria and the delay-independent stability criteria. In many literatures, linear matrix inequality(LMI) method is often used to obtain the delay-dependent stability criteria for GRNs with time delay, which has very complex computation. In order to reduce complex computation, the Dini derivative method is applied in [20] to obtain a stability criterion for a class of the continuous-time GRNs. This obtained stability criterion is given in the form of several elementary inequalities, which are very simple to be easily implemented on computer. The main purpose of this paper is to improve the Dini derivative method proposed in [20] by constructing a new Lyapunov functional. In this paper, a new Lyapunov functional is constructed to obtain a sufficient condition of the global exponential stability for a class of the continuous-time GRNs with time delay, which is given in the form of several elementary inequalities. For the continuous-time GRNs with delays, the semi-discretization technique and the IMEX-θ method are used to obtain two discretetime GRN systems, which can preserve the dynamical characteristics of the continuous-time systems. Furthermore, it is shown that under the same sufficient conditions obtained earlier, these two discrete-time GRN systems are also globally exponentially stable. Comparing with [20, Assumption 1], the global exponential stability criterion obtained in this paper is less conservative, which has been shown in Remark 2 and Example 2. The remaining part of the paper is organized as follows. In Section 2.1, the model of continuous-time GRNs is introduced. Based on this model, a sufficient condition of the global exponential stability for continuous-time GRNs is investigated in Section 2.2. By using the semi-discretization technique and the IMEX-θ method, two kinds of discrete-time systems corresponding to the continuous-time GRN model are derived in Sections 3.1 and 3.2, respectively. It is shown that under the same conditions, the exponential convergence characteristics are also preserved by these two discrete-time analogs. In Section 3.3, two numerical examples are given to demonstrate the effectiveness of the obtained theoretical results. Finally, some conclusions are drawn in Section 4. 2. Global exponential stability of continuous-time GRNs 2.1. Model description and preliminaries In this paper, R and IN denote the 1-dimensional Euclidean space and the index set {1, 2 . . . , N}, respectively. C ([a, b], RN ) represents the set of all N-dimensional continuous functions defined on the interval [a, b], and [a, b]N means [a, b] × . . . × [a, b].







N

This paper considers the following delayed GRNs in [20]:



m˙ i (t ) = −ai mi (t ) +

N j=1

bi j f j ( p j (t − τi j )) + Ii ,

p˙ i (t ) = −ci pi (t ) + di mi (t − σi ),

(1)

i ∈ IN ,

where mi (t) and pi (t) are the concentrations of the mRNA and the protein of the ith node at time t, respectively; ai > 0 and ci > 0 represent, respectively, the degradation rates of the mRNA and the protein; di > 0 is the translation rate of the ith gene; τ ij > 0 and σ i > 0 denote the feedback regulation delay and the translation delay, respectively; bij is defined as follows: (i) If there is no connection between j and i, then bi j = 0, (ii) If transcription factor j activates gene i, then bi j = αi j , (iii) If transcription factor j  represses gene i, then bi j = −αi j ;Ii = j=S αi j is the basal rates, in which Si is the set of repressors of the ith gene. The nonlinear i

feedback regulation function fj (pj (t)) possesses the following form: f j ( p j (t )) = ( p j (t )/β j )h j /(1 + ( p j (t )/β j )h j ), where hj is the Hill coefficient, and β j > 0 is a constant. Obviously, fj (·) ( j ∈ IN ) are monotonically nondecreasing functions with saturation, it satisfies

f j (0) = 0,

0≤

f j (ν1 ) − f j (ν2 ) ≤ hj ν1 − ν2

(2)

for all ν1 , ν2 ∈ R with ν 1 = ν 2 , and | f j (ν1 )| ≤ N j < +∞, j ∈ IN , where hj > 0 and Nj > 0 are constants; The initial condition of GRNs (1) is given as follows: mi (s) = φi (s), s ∈ [−σ , 0], pi (s) = ψi (s), s ∈ [−τ , 0], where φi ( · ) ∈ C ([−σ , 0], R), ψi ( · ) ∈ C ([−τ , 0], R), τ = maxi, j∈IN {τi j } and σ = maxi∈IN {σi }.  For = maxi∈IN {di ( Nj=1 |bi j |N j + |Ii |)/(ai ci )}, define a continuous mapping T: [0, ϖ]N → [0, ϖ]N as follows: T ( p) =  (T1 ( p), T2 ( p), . . . , TN ( p))T with Ti ( p) = adci Nj=1 bi j f j ( p j ) + adci Ii , i ∈ IN , for any p ∈ [0, ϖ]N . According to Brouwer’s fixed i i

i i

point theorem [21,22], there exists one fixed point p∗ for the continuous mapping T. Let m∗ = (m∗1 , m∗2 , . . . , m∗N )T and p∗ = ( p∗1 , p∗2 , . . . , p∗N )T . Then



 −ai m∗i + Nj=1 bi j f j ( p∗j ) + Ii = 0, ∗ −ci pi + di m∗i = 0, i ∈ IN ,

and thus (m∗ , p∗ ) is the unique equilibrium point of GRNs (1). Definition 1. [23] For given a function V (t ) : [t0 , +∞) → R, the upper right Dini derivative of V(t), D+V (t ), is

D+V (t ) = lim sup h→0+

V (t + h) − V (t ) . h

(3)

630

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

Remark 1. Similar to the definition of the upper right Dini derivative, one can define the lower right Dini derivative, the upper left Dini derivative and the lower left Dini derivative. These four kinds of derivatives are unify called Dini derivative. The Dini derivative is a generalization of derivative, i.e., when the derivative of a function exists, its Dini derivative must exist and is equal to the derivative, but the opposite may not be true. The Dini derivative can be applied to the continuous and non-differentiable functions, such as Lyapunov functional with absolute value symbol. Thus, Dini derivative of V(t) in (10) can be used to analyze the stability of GRN (1), but its derivative does not exist (see [24,25]) .  Definition 2. GRNs (1) is said to be globally exponentially stable if there exist constants α > 0 and ξ > 0 such that N i=1 {|mi (t ) − N ∗ ∗ − α t ∗ for all t ≥ 0, where Ki = sups∈[−σ ,0] |φi (s) − mi | and Li = sups∈[−τ ,0] |ψi (s) − p∗i |. mi | + | pi (t ) − pi |} ≤ ξ i=1 {Ki + Li }e 2.2. Global exponential stability of continuous-time GRNs In this section, we will investigate the globally exponential stability of continuous-time GRNs (1) with delays, a sufficient condition is obtained immediately. Theorem 1. If

ai ci > di hi

N 

|bi j |, ∀i ∈ IN ,

(4)

j=1

then GRNs (1) is globally exponentially stable.  Proof. From (4), there exists a constant γ i > 0 such that γi ai − di > 0, ci − γi hi Nj=1 |b ji | > 0 for every i ∈ IN . Define the contin uous functions: Fi (μ) = γi ai − γi μ − di eμσi , Gi (ν) = ci − ν − γi hi Nj=1 |b ji |eντ ji . It is easy to verify that

Fi (0) = γi ai − di > 0,

Gi (0) = ci − γi hi

N 

|b ji | > 0,

j=1

and

limμ→+∞ Fi (μ) = −∞,

limν →+∞ Gi (ν) = −∞.

According to the properties of the continuous functions, there exist constants μi , νi ∈ (0, +∞) such that Fi (μi ) = γi ai − γi μi −  di eμi σi = 0,Gi (νi ) = ci − νi − γi hi Nj=1 |b ji |eνi τ ji = 0. Furthermore, we have

F˙i (μ) = −γi − di σi eμσi < 0 and

G˙ i (ν) = −1 − γi hi

N 

τ ji |b ji |eντ ji < 0.

j=1

This means that Fi (μ) and Gi (ν ) are strictly monotone decreasing functions. Thus,

Fi (μ) ≥ 0,

0 ≤ μ ≤ μi ; Fi (μ) < 0,

μi < μ < +∞

and

Gi (ν) ≥ 0,

0 ≤ ν ≤ νi ; Gi (ν) < 0,

νi < ν < +∞.

Let

α = min {μ1 , μ2 , . . . , μN , ν1 , ν2 , . . . , νN }. Then

{Fi (μ) ≥ 0, μ ∈ [0, α ], ∀i ∈ IN , Gi (ν) ≥ 0, ν ∈ [0, α ], ∀i ∈ IN .

(5)

Suppose (mi (t), pi (t)) and (m∗i , p∗i ) are an arbitrary solution and the equilibrium point of GRNs (1), respectively. Define

xi (t ) = mi (t ) − m∗i ,

yi (t ) = pi (t ) − p∗i

and

f˜j (y j (t − τi j )) = f j (y j (t − τi j ) + p∗j ) − f j ( p∗j ). Then, GRNs (1) can be rewritten as:

{x˙ i (t ) = −ai xi (t ) +

N  j=1

bi j f˜j (y j (t − τi j )), y˙ i (t ) = −ci yi (t ) + di xi (t − σi ),

i ∈ IN .

(6)

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

631

Set

Mi (t ) = eαt xi (t ),

t ∈ [−σ , +∞)

and

Pi (t ) = eαt yi (t ),

t ∈ [−τ , +∞).

Applying the definition of the upper right Dini derivative (3), we will discuss the upper right Dini derivatives of |Mi (t)| and |Pi (t)|, respectively. The following two cases are considered. Case 1: Mi (t) = 0. Using (6), we have +

d |Mi (t )| = sgn(xi (t ))(α eαt xi (t ) + eαt x˙ i (t )) dt N  ≤ (α − ai )|Mi (t )| + |bi j |h j eατi j |Pj (t − τi j )|.

(7)

j=1

Case 2: Mi (t ) = 0. By using (6), we obtain +

d |Mi (t )| = dt ≤

α eαt |xi (t )| + eαt |x˙ i (t )| (α − ai )|Mi (t )| +

N 

|bi j |h j eατi j |Pj (t − τi j )|.

(8)

j=1

Combining the two cases above, we can conclude

 d |Mi (t )| ≤ (α − ai )|Mi (t )| + |bi j |h j eατi j |Pj (t − τi j )|. dt +

N

(9a)

j=1

Similar to (9a), it is easy to deduce from (6) that +

d |Pi (t )| ≤ (α − ci )|Pi (t )| + di eασi |Mi (t − σi )|. dt

(9b)

Construct a Lyapunov functional candidate as follows:

V (t ) =

N 

(V1i (t ) + V2i (t )), t ≥ 0.

(10)

i=1

with

V1i (t ) = γi |Mi (t )| + γi

N 

|bi j |h j eατi j

j=1

and

V2i (t ) = |Pi (t )| + di eασi



t

t−σi



t

t−τi j

|Pj (s)|ds

|Mi (s)|ds.

Combining (5) with (9), we obtain

 d |V (t )| ≤− (Fi (α)|Mi (t )| + Gi (α)|Pi (t )|) ≤ 0. dt +

N

(11)

i=1

Let γ = min{γ1 , γ2 , . . . , γN , 1}, then γ > 0. In view of (10) and (11), it is clear that

γ

N 

{|Mi (t )| + |Pi (t )|} ≤ V (t ) ≤ V (0).

(12)

i=1

Let

Ki = sup

s∈[−σ ,0]

and

|xi (s)|, Li = sup |yi (s)|, ξ0 = max{ξi }, s∈[−τ ,0]

ξi = max γi + di eασi σi , 1 + γi hi

i∈IN

N  j=1

 |b ji |eατ ji τ ji .

632

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

Then

V (0) ≤ ξ0

N 

(Ki + Li ).

(13)

i=1

This, together with (12), implies that N 

(|mi (t ) − m∗i | + | pi (t ) − p∗i |) ≤

i=1

N ξ0  (K + L )e−αt , γ i=1 i i

From Definition 2, GRNs (1) is globally exponentially stable. This proof is completed.



Remark 2. If [20, Assumption 1] is satisfied, then it is obvious that (4) holds. However, the opposite is not true (see Example 2 below). This implies the sufficient condition given in Theorem 1 is weaker than [20, Theorem 1]. 3. Global exponential stability of discrete-time GRNs 3.1. Semi-discretization model In this section, the semi-discretization method is applied to the ordinary differential equations. It has been proved in [15] that this kind of method can preserve the dynamical characteristics of the continuous-time systems to some extent. Suppose the discretization step-size h is a fixed positive real number, and for any t ∈ [0, +∞), there exists a non-negative integer n such that t ∈ [nh, (n + 1)h), where n = [t/h], and [ · ] represents the integral function. Let [τi j /h] = ωi j , [σi /h] = κi , i, j ∈ IN . Then, the system (1) can be transformed into the following approximated form:



 m˙ i (t ) = −ai mi (t ) + Nj=1 bi j f j ( p j (n − ωi j )) + Ii p˙ i (t ) = −ci pi (t ) + di mi (n − κi ),

(14)

where p j (n − ωi j ) and mi (n − κi ) represent pi (nh − ωi j h) and mi (nh − κi h), respectively. Let

θi (h) = (1 − e−ai h )/ai > 0, ϕi (h) = (1 − e−ci h )/ci > 0 Multiplying both sides of (14) with eai t , integrating over the interval [nh, t] with t < (n + 1)h, and taking the limit t → (n + 1)h, we obtain



N

mi (n + 1) = e−ai h mi (n) + θi (h) j=1 bi j f j ( p j (n − ωi j )) + Ii pi (n + 1) = e−ci h pi (n) + ϕi (h)(di mi (n − κi )).

(15)

where mi (n) and pi (n) are used to represent mi (nh) and pi (nh), respectively. The initial condition of system (15) is given as

{mi (l ) = φi (l ), l ∈ [−κ , 0], pi (l ) = ψi (l ), l ∈ [−ω, 0], where φ i (l) and ψ i (l) stand for φ i (lh) and ψ i (lh), which are all bounded functions. After the simple calculation, the discrete-time system (15) and the continuous-time system (1) have the same equilibrium. It has been shown when h → 0, the discrete-time model (15) can be reduced to the original continuous-time system (1). Definition 3. The GRNs (15) is said to be globally exponentially stable if there exist constants λ > 1 and δ > 0 such that N 

{|mi (n) − m∗i | + | pi (n) − p∗i |} ≤ δ

i=1

holds for all n ∈

N 

{Ki + Li }λ−n

i=1

Z+ ,

where Ki = supl∈[−κ ,0]Z |φi (l ) − m∗i | and Li = supl∈[−ω,0]Z |ψi (l ) − p∗i |.

Theorem 2. If the condition (4) is satisfied, then the system (15) has a unique equilibrium (m∗T , p∗T ) which is globally exponentially stable. Proof. Since the discrete-time GRN model (15) has the same equilibrium with the continuous-time system (1), it follows from the discussions in Section 2.1 and Theorem 1 that the equilibrium of system (15) does exist and is unique. In the following, we only need to prove that the equilibrium (m∗T , p∗T ) is globally exponentially stable.  From (4), there exists a constant γ i > 0 such that γi ai − di > 0, ci − γi hi Nj=1 |b ji | > 0. Define the continuous functions:

Fi (ςi ) = γi − γi ςi e−ai h − di θi (h)ςiκi +1 and

Gi (ζi ) = 1 − ζi e−ci h − γi hi ϕi (h)

N  j=1

|b ji |ζi ω ji +1 .

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

N

It is easy to verify Fi (1) = θi (h)[γi ai − di ] > 0, Gi (1) = ϕi (h)[ci − γi hi

lim Fi (ςi ) = −∞,

ςi →+∞

j=1

633

|b ji |] > 0, and

lim Gi (ζi ) = −∞.

ζi →+∞

∗

∗

According to the properties of the continuous functions, there must exist ς i ∗ and ζi ∈ (1, +∞) such that Fi (ςi ∗ ) = Gi (ζi ) = 0, Fi (ςi ) > 0, ςi ∈ (1, ςi ∗ ) and Gi (ζi ) > 0, ζi ∈ (1, ζi∗ ). Let

λ = min{ς1∗ , ς2∗ , . . . , ςN∗ , ζ1∗ , ζ2∗ , . . . , ζN∗ } It is clear λ > 1. In the following, define a Lyapunov functional candidate as

V2 (n) =

N 



γi Mi (n) + γi

n−1 N  

|bi j |h j ϕ j (h)λωi j +1 Pj (l ) + Pi (n) +

j=1 l=n−ωi j

i=1

n−1 

 di θi (h)λκi +1 Mi (l )

(16)

l=n−κi

where Mi (n), Pi (n) are defined as

Mi (n) = λn |mi (n) − mi ∗ |θi−1 (h),

n ∈ [−κ , +∞)Z ,

(17a)

and

Pi (n) = λn | pi (n) − pi ∗ |ϕi−1 (h),

n ∈ [−ω, +∞)Z .

(17b)

From the GRN model (15), we can get

|mi (n + 1) − mi ∗ | ≤ e−ai h |mi (n) − mi ∗ | + θi (h)

N 

|bi j ||h j || p j (n − ωi j ) − p j ∗ |

(18)

j=1

Similar to (18), it is easy to obtain

| pi (n + 1) − pi ∗ | ≤ e−ci h | pi (n) − pi ∗ | + ϕi (h)di |mi (n − κi ) − mi ∗ |

(19)

From (17), (18) and (19), it follows

Mi (n + 1) ≤ λe−ai h Mi (n) +

N 

|bi j |h j ϕi (h)λωi j +1 Pj (n − ωi j ),

j=1

and

Pi (n + 1) ≤ λe−ci h Pi (n) + di θi (h)λκi +1 Mi (n − κi ). where i ∈ IN , n ∈ Z0 + . In the following, we shall compute the difference along the solutions of system (15).

V2 (n) = V2 (n + 1) − V2 (n)  ≤−

N 

(γi − γi λe

−ai h

− di θi (h)λκi +1 )Mi (n) + (1 − λe−ci h − γi

i=1

=−

N 

N 

 |b ji |hi ϕi

(h)λω ji +1 )P (n)

j=1

(Fi (λ)Mi (n) + Gi (λ)Pi (n))

i=1

≤0 This, together with (16) yields

λn γ

N 

N 

i=1

i=1

{|mi (n) − mi ∗ |θi−1 (h) + | pi (n) − pi ∗ |ϕi−1 (h) = γ

where γ = min{γ1 , γ2 , . . . , γN , 1} > 0 and

V2 (0) ≤

N  i=1

Thus



1 + di κi λκi +1 )Ki + (γi θi (h)



{Mi (n) + Pi (n)} ≤ V2 (n) ≤ V2 (0),  

 1 + γi |b ji |hi λω ji +1 ω ji Li ϕi (h) N

j=1

i

634

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643 N 

{|mi (n) − m∗i | + | pi (n) − p∗i |} ≤

i=1

N δ0  {K + Li }λ−n , γ i=1 i

where

δ0 =

max{maxi∈IN {θi (h)}, maxi∈IN {ϕi (h)}} max{δ1 , δ2 } min{mini∈IN {θi (h)}, mini∈IN {ϕi (h)}}

δ1 = max{γi + di θi (h)κi λκi +1 } ≥ γi i∈IN

and

 δ2 = max 1 + γi ϕi (h) i∈IN

N 

 |b ji |hi λω ji +1 ω ji

≥ 1.

j=1

From Definition 3, it follows that system (15) is globally exponentially stable. This proof is completed.



3.2. Implicit-explicit-θ discretization model In this subsection, the implicit-explicit-θ discretization method is used to obtain the following discrete-time GRN model:





N  1 − (1 − θ )hai h {mi (n + 1) = mi (n) + bi j f j ( p j (n − ωi j )) + Ii , 1 + θ hai 1 + θ hai j=1

1 − (1 − θ )hci di h pi (n) + mi (n − κi ), pi (n + 1) = 1 + θ hci 1 + θ hci

θ ∈ [0, 1].

(20)

where h, ωij and κ i are the same as defined in (15). The initial condition of system (20) is given as one of system (15). The equilibrium of system (20) is denoted as (m∗T , p∗T ), where m∗ = (m∗1 , m∗2 , . . . , m∗N )T and p∗ = ( p∗1 , p∗2 , . . . , p∗N )T . After the simple computation, it follows that the discrete-time system (20) and the continuous-time system (1) have the same equilibrium (m∗T , p∗T ). Theorem 3. Under the hypothesis of Theorem 1, the system (20) has a unique equilibrium (m∗T , p∗T ), which is globally exponentially stable. Proof. Similar to the proof of Theorem 2, this proof can be obtained immediately.  Remark 3. In many literatures, researchers have studied the problems of model by the semi-discretization method [26–28]; In different kinds of discrete-time models, the IMEX-θ discretization method is widely employed [29,30]. By using the two kinds of methods, the continuous-time systems can be transformed into two discrete-time GRN systems, respectively. The two discretetime GRNs obtained have two advantages: (i) They can preserve the dynamical characteristics of the continuous-time systems; (ii) They are easy to implement on the computer. In addition, the convergence rate of the discrete-time GRN derived from the IMEX-θ method can be improved by choosing the appropriate parameter. Comparing with LMI method, the semi-discretization technique and the IMEX-θ method have less computation, which is simple to be easily implemented on computer. 3.3. Illustrative examples In this section, we provide two numerical examples to demonstrate the effectiveness of the theoretical results obtained. Example 1. Consider GRNs (1) and (15) with

a1 = 5.1,

a2 = 4.6,

a3 = 5.5,

c1 = 7.1,

c2 = 7.8,

c3 = 7.3,

d1 = 0.8,

d2 = 0.9,

d3 = 0.7,

b11 = b12 = b22 = b23 = b31 = b33 = 0, b13 = b21 = b32 = −2.5, I1 = I2 = I3 = 2.5, f j ( p j ) = 0.7p j 2 /(1 + p j 2 ), h j = 0.65 ( j = 1, 2, 3). This example is taken from [31]. After the simple computation, it is found that GRNs (1) and (15) satisfy the condition (4) have the unique equilibrium point (m∗ , p∗ )T , where

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

635

0.8 m1(t) m2(t) m3(t)

0.75

mRNA concentration

0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3

0

1

2

3

4

5

Times Fig. 1. Trajectories of mRNA concentration mi (t).

0.8 p1(t) p2(t) p3(t)

0.6

protein concentration

0.4

0.2

0

−0.2

−0.4

−0.6

0

1

2

3

4

5

Times Fig. 2. Trajectories of protein concentration pi (t).

m∗ = (0.4895, 0.5423, 0.4523)T , p∗ = (0.0552, 0.0626, 0.0435)T . From Theorems 1 and 2, it follows that the considered GRNs (1) and (15) are globally exponentially stable. To further illustrate the theoretical result, choose the initial conditions φ1 ≡ (0.3, 0.8, 0.5)T , ψ1 ≡ (0.2, 0.7, −0.5)T and the time delays τi j = 0.3, σi = 0.4 (i, j ∈ I3 ). By using MATLAB, the trajectories of the continuous-time GRN (1) are shown in Figs. 1 and 2, and the trajectories of the discrete-time system (15) are shown in Figs. 3 and 4. From Figs. 1–4, we can observe the mRNA and protein concentrations (i.e., the system states) are stable.

636

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

0.8 m1(t) m2(t) m3(t)

0.75

mRNA concentration

0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3

0

2

4

6

8

10

Times Fig. 3. Trajectories of mRNA concentration mi (t).

0.8 p1(t) p2(t) p3(t)

0.6

protein concentration

0.4

0.2

0

−0.2

−0.4

−0.6

0

2

4

6

8

10

Times Fig. 4. Trajectories of protein concentration pi (t).

In Theorem 1, the condition (4) is delay-independent, this means that the continuous-time GRN (1) and the discrete-time system (15) are both globally exponentially stable for arbitrary time delay when the condition (4) is satisfied. To give further explanation for large delay case, we choose τi j = σi = 100 (i, j ∈ I3 ) and other data above. By using MATLAB, the trajectories of the continuous-time GRN (1) are shown in Figs. 5 and 6, and the trajectories of the discrete-time system (15) are shown in Figs. 7 and 8. From these figures, we can observe the mRNA and protein concentrations (i.e., the system states) are stable. Compared with Figs. 1−4, we can see that, in Figs. 5–8, the larger value the time delay, the longer time the GRNs reach steady state. However, all trajectories in Figs. 1–8 are stable.

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

637

0.8 m1(t) m2(t) m3(t)

0.75

mRNA concentration

0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3

0

500

1000

1500

Times Fig. 5. Trajectories of protein concentration mi (t).

0.5 p1(t) p2(t) p3(t)

0.45

protein concentration

0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

0

500

1000 Times

Fig. 6. Trajectories of protein concentration pi (t).

Example 2. Consider GRNs (1) and (15) with

a1 = 0.28,

a2 = 0.32,

a3 = 0.40,

c1 = 0.45,

c2 = 0.42,

c3 = 0.35,

d1 = 0.30,

d2 = 0.35,

d3 = 0.42,

b11 = b12 = b22 = b23 = b31 = b33 = 0, b13 = b21 = b32 = −0.45,

1500

638

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

0.8 m1(t) m2(t) m3(t)

0.75

mRNA concentration

0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3

0

500

1000 Times

1500

2000

Fig. 7. Trajectories of protein concentration mi (t).

0.5 p1(t) p2(t) p3(t)

0.45

protein concentration

0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

0

500

1000 Times

1500

2000

Fig. 8. Trajectories of protein concentration pi (t).

I1 = I2 = I3 = 0.45, f j ( p j ) = p j 2 /(1 + p j 2 ), h j = 0.65( j = 1, 2, 3). These data are taken from [20] except ai , ci and di , i ∈ IN . After simple computation, it is found that GRNs (1) and (15) have the unique equilibrium point (m∗ , p∗ )T , where

m∗ = (0.9757, 0.9882, 0.6704)T , p∗ = (0.6505, 0.8235, 0.8045)T .

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

639

1.3 m1(t) m2(t) m3(t)

1.2

mRNA concentration

1.1 1 0.9 0.8 0.7 0.6 0.5

0

0.5

1 Times

1.5

2 4

x 10

Fig. 9. Trajectories of protein concentration mi (t).

2 p1(t) p2(t) p3(t)

1.8

protein concentration

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

0.5

1 Times

1.5

2 4

x 10

Fig. 10. Trajectories of protein concentration pi (t).

It can easily be verified that [20, Assumption 1] is not true, and hence one can not test the global exponential stability of the considered GRNs by using [20, Theorem 1–3]. However, the condition (4) of this paper is satisfied. It follows from Theorems 1 and 2 that the GRNs (1) and (15) are globally exponentially stable. To simulate the theoretical results, choose the initial conditions φ 2 ≡ (0.62, 0.52, 0.88)T , ψ 2 ≡ (0.81, 0.43, 0.67)T and the time delays τi j = σi = 100 (i, j ∈ I3 ) and using MATLAB, the trajectories of the continuous-time GRN (1) are shown in Figs. 9 and 10, and the trajectories of the discrete-time system (15) are shown in Figs. 11 and 12. From Figs. 9–12, we can observe the mRNA and protein concentrations (i.e., the system states) are stable. In order to further show large delay case, we choose the time-delay τi j = σi = 300 (i, j ∈ I3 ), by using MATLAB, the trajectories of the continuous-time GRN (1) are shown in Figs. 13 and 14, and the trajectories of the discrete-time system (15) are shown

640

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

1.3 m1(t) m2(t) m3(t)

1.2

mRNA concentration

1.1 1 0.9 0.8 0.7 0.6 0.5

0

0.5

1 Times

1.5

2 4

x 10

Fig. 11. Trajectories of protein concentration mi (t).

2 p1(t) p2(t) p3(t)

1.8

protein concentration

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

0.5

1 Times

1.5

2 4

x 10

Fig. 12. Trajectories of protein concentration pi (t).

in Figs. 15 and 16. From these figures, we can observe the mRNA and protein concentrations (i.e., the system states) are stable. Compared with Figs. 9–12, we can see the larger value the time delay, the longer time the GRNs reach steady state. 4. Conclusions In this paper, Dini derivative method is used to analyze the global exponential stability of the continuous-time genetic regulatory networks (GRNs) with delays. Based on the Lyapunov stability theory and some inequality techniques, the same sufficient condition is derived to guarantee the globally exponentially stability for both the continuous-time GRN systems and the discrete-time GRN systems, where the discrete-time GRNs corresponding to the continuous-time GRNs have been

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

641

1.3 m1(t) m2(t) m3(t)

1.2

mRNA concentration

1.1 1 0.9 0.8 0.7 0.6 0.5

0

0.5

1 Times

1.5

2 4

x 10

Fig. 13. Trajectories of protein concentration mi (t).

2 p1(t) p2(t) p3(t)

1.8

protein concentration

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

0.5

1 Times

1.5

2 4

x 10

Fig. 14. Trajectories of protein concentration pi (t).

derived by the semi-discretization technique and the IMEX-θ method. This sufficient condition is less conservative than one obtained in [20]. Two numerical examples are provided to show the effectiveness of the proposed stability analysis approach. It is worth mentioning that the stability criteria obtained in this paper has a certain conservativeness, since it is delayindependent. This motivates us to improve the Dini derivative method proposed in this paper to obtain delay-dependent stability criteria of the delayed GRNs in the near future. In addition, there are still several related problems to need be solved, such as the robust stability problem (see [32]) and the state estimation problem (see [33–35]).

642

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

1.3 m1(t) m2(t) m3(t)

1.2

mRNA concentration

1.1 1 0.9 0.8 0.7 0.6 0.5

0

0.5

1 Times

1.5

2 4

x 10

Fig. 15. Trajectories of protein concentration mi (t).

2 p1(t) p2(t) p3(t)

1.8

protein concentration

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

0.5

1 Times

1.5

2 4

x 10

Fig. 16. Trajectories of protein concentration pi (t).

Acknowledgments The authors thank the editor and the referees for their valuable suggestions to improve the quality of this paper. This work was partially supported by the National Natural Science Foundation of China (no. 11371006), Natural Science Foundation of Heilongjiang Province (no. A201416), the Undergraduate Innovation and Entrepreneurship Program of Heilongjiang University (no. 2014SX04) and the Laboratory Project of Heilongjiang University. References [1] W. Zhang, J. Fang, Y. Tang, Robust stability for genetic regulatory networks with linear fractional uncertainties, Commun. Nonlinear Sci. Numer. Simul. 17 (4) (2012) 1753–1765.

T. Liu et al. / Applied Mathematics and Computation 274 (2016) 628–643

643

[2] W. Wang, S. Zhong, Delay-dependent stability criteria for genetic regulatory networks with time-varying delays and nonlinear disturbance, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 3597–3611. [3] H.D. Jong, Modeling and simulation of genetic regulatory systems: a literature review, J. Comput. Biol. 9 (1) (2002) 67–103. [4] J.L. Liu, D. Yue, Asymptotic and robust stability of T-S fuzzy genetic regulatory networks with time-varying delays, Int. J. Robust Nonlinear Control 22 (8) (2012) 827–840. [5] N. Friedman, M. Linial, I. Nachman, Using Bayesian networks to analyze expression data, J. Comput. Biol. 7 (2000) 601–620. [6] S. Hardy, P.N. Robillard, Modeling and simulation of molecular biology systems using Petri nets: modeling goals of various approaches, J. Bioinform. Comput. Biol. 2 (04) (2004) 619–637. [7] R. Somogyi, C. Sniegoski, Modeling the complexity of genetic networks: understanding multigenic and pleiotropic regulation, Complexity 1 (1996) 45–63. [8] P. Smolen, D.A. Baxter, J.H. Byrne, Mathematical modeling of gene networks, Neuron 26 (3) (2000) 567–580. [9] L. Chen, K. Aihara, Stability of genetic regulatory networks with time delay, IEEE Trans. Circuits Syst. Regul. Pap. 49 (5) (2002) 602–608. [10] P. Smolen, D.A. Baxter, J.H. Byrne, Modeling circadian oscillations with interlocking positive and negative feedback loops, J. Neurosci. 21 (17) (2001) 6644– 6656. [11] Y.T. Wang, X. Zhang, Z.R. Hu, Delay-dependent robust H∞ filtering of uncertain stochastic genetic regulatory networks with mixed time-varying delays, Neurocomputing 166 (2015) 346–356. [12] L.L. Chen, Y. Zhou, X. Zhang, Guaranteed cost control for uncertain genetic regulatory networks with interval time-varying delays, Neurocomputing 131 (2014) 105–112. [13] N.A.M. Monk, Oscillatory expression of Hes1, p53, and NF-κ B driven by transcriptional time delays, Curr. Biol. 13 (16) (2003) 1409–1413. [14] D. Zhang, L. Yu, Passivity analysis for stochastic Markovian switching genetic regulatory networks with time-varying delays, Commun. Nonlinear Sci. Numer. Simul. 16 (8) (2011) 2985–2992. [15] M. Xiao, J.D. Cao, Genetic oscillation deduced from Hopf bifurcation in a genetic regulatory network with delays, Math. Biosci. 215 (1) (2008) 55–63. [16] X. Zhang, A.H. Yu, G.D. Zhang, M-matrix-based delay-range-dependent global asymptotical stability criterion for genetic regulatory networks with timevarying delays, Neurocomputing 113 (2013) 8–15. [17] X. Zhang, L.G. Wu, S.C. Cui, An improved integral to stability analysis of genetic regulatory networks with interval time-varying delays, IEEE/ACM Trans. Comput. Biol. Bioinform. 19 (13) (2015) 1493–1518. [18] Y.T. Wang, A.H. Yu, X. Zhang, Robust stability of stochastic genetic regulatory networks with time-varying delays: a delay fractioning approach, Neural Comput. Appl. 23 (5) (2013) 1217–1227. [19] Y.Y. Han, X. Zhang, Y.T. Wang, Asymptotic stability criteria for genetic regulatory networks with time-varying delays and reaction–diffusion terms, Circuits Syst. Signal Process. 34 (2015) 3161–3190. [20] J.Q. Hu, J.L. Liang, J.D. Cao, Stability analysis for genetic regulatory networks with delays: the continuous-time case and the discrete-time case, Appl. Math. Comput. 220 (2013) 507–517. [21] M.M. Day, Fixed-point theorems for compact convex sets, Ill. J. Math. 5 (1961) 585–590. [22] L.E.J. Brouwer, Über abbildung yon mannigfaltigkeiten, Math. Ann. 71 (1) (1911) 97–115. [23] V.B. Kolmanovskii, A. Myshkis, Applied Theory of Functional Differential Equations, Kluwer Academic Publishers, Boston, 1992. [24] C. Briat, H. Hjalmarsson, K.H. Johansson, U.T. Jonsson, G. Karlsson, H. Sandberg, Nonlinear state-dependent delay modeling and stability analysis of internet congestion control, in: Decision and Control (CDC), 2010, 49th IEEE Conference on, 2010, pp. 1484–1491. [25] M. Wang, T. Zhou, H. Fang, Global exponential stability of mam neural network with time-varying delays, in: Computational Intelligence and Software Engineering (CISE), 2010 International Conference on, 2010, pp. 1–4. [26] T. Insperger, G. Stpn, Updated semi-discretization method for periodic delay-differential equations with discrete delay, Int. J. Numer. Methods Eng. 61 (1) (2004) 117–141. [27] A. Chanda, A. Fischer, Stability analysis of a thin-walled cylinder in turning operation using the semi-discretization method, Acta Mechanica Sinica 30 (2) (2014) 214–222. [28] O. Elbeyli, J.Q. Sun, Efficiency and accuracy of semi-discretization schemes for time delayed feedback control systems, in: ASME 2003 International Mechanical Engineering Congress and Exposition, 2003, pp. 1019–1036. [29] T. Koto, Stability of imex rungeckutta methods for delay differential equations, J. Comput. Appl. Math. 211 (2008) 201–212. [30] W.S. Wang, S.F. Li, Nonlinear stability of rungeckutta methods for neutral delay differential equations, J. Comput. Appl. Math. 214 (1) (2008) 175–185. [31] X. Lou, Q. Ye, B. Cui, Exponential stability of genetic regulatory networks with random delays, Neurocomputing 73 (2010) 759–769. [32] S. Lakshmanan, F.A. Rihan, Stability analysis of the differential genetic regulatory networks model with time-varying delays and markovian jumping parameters, Nonlinear Anal.: Hybrid Syst. 14 (2014) 1–15. [33] T.H. Lee, S. Lakshmanan, J.H. Park, P. Balasubramaniam, State estimation for genetic regulatory networks with mode-dependent leakage delays, time-varying delays, and Markovian jumping parameters, NanoBiosci. IEEE Trans. 12 (4) (2013) 363–375. [34] S. Lakshmanan, H. Park, H.Y. Jung, P. Balasubramaniam, S.M. Lee, Design of state estimator for genetic regulatory networks with time-varying delays and randomly occurring uncertainties, Biosystems 111 (1) (2013) 51–70. [35] R. Sakthivel, K. Mathiyalagan, S. Lakshmanan, J.H. Park, Robust state estimation for discrete-time genetic regulatory networks with randomly occurring uncertainties, Nonlinear Dyn. 74 (4) (2013) 1297–1315.