Stability analysis around singular stationary points and computer algebra systems Текст научной статьи по специальности «Математика»

YJ\K 517.929.4, 51-76

STABILITY ANALYSIS AROUND SINGULAR STATIONARY POINTS AND COMPUTER ALGEBRA SYSTEMS © Irina Shlykova, Arcady Ponosov, Yuriy Nepomnyashchikh and Andrei Shindiapin

Key words: gene regulatory systems, step functions, stability of non-continuous systems.

A CAS-based method to study asymptotic properties of solutions to systems of differential equations with distributed time-delays and Boolean-type nonlinearities (step functions) is offered. Such systems arise in many applications, but this paper deals with specific examples of such systems coming from genetic regulatory networks. A challenge is to analyze stable stationary points which belong to the discontinuity set of the system (thresholds). In this work we do the stability analysis around such points in the presence of delays. The basic technical tool consists in replacing step functions by the so-called "logoid functions" and investigating the smooth systems thus obtained.

1. Introduction

We study asymptotically stable steady states (stationary points) of the system

Xi = Fi(Z\,Zn) Zn)xi,

Zi — Zi{yi),

W(i) = («*!)(*) (t> 0),

Vi = Xi (i = 2, ...,n).

This system describes a specific gene regulatory network with autoregulation.

The functions Fi, Gi are affine in each Zi and satisfy

Fi(Z\,..., Zn) > 0, Gi(Z\,Zn) > 0

(0 < Zi < 1, i = 1, ...,n). The input variables yt are described by nonlinear Volterra ("delay") operators.

Below we assume that 5? is the integral operator given by

(Ux){t) = c0x{t)+ f K(t-s)x(s)ds, 0, (2)

J—OO

where K(u) = cvKv(u), c„ ^ 0 (v = 0,1, Q) + Ylv=i Cv ~ ^ an(^

(3)

Each Zi is a steep sigmoid function depending on the input variable yi, i.e. Z^yi) — T,(yi,9i,qi). In the vicinity of the threshold value Qi the response function Zt switches from 0 to 1. Thus, in the limit the response function is close to the step function having the unit jump at i^ = &i.

The aim is to examine stability properties of SSP when a critical variable is subject to a delay. The main method we use is a sort of localization procedure for SSP presented in [6] for rather general genetic regulatory networks with delay.

Вестник ТГУ, т. 15, вып.2,2010

2. Response functions ~~

Let describe the properties of general logoid functions. These will serve as response functions in the model given by (1).

Let Z = E(y, 9, q) be any function defined for у Є R, 9 > 0, 0 < q < q° .

Assumption 1: T,(y,9,q) is continuous in (y,q) Є R x (0, q°) for all в > 0, continuously differentiable w.r.t. у Є R for all 9 > 0,0 < q < q°, and ^H(y,9,q) > 0 on the set {у Є R : 0 < Z(y,6,q) < 1} .

Assumption 2: E(y, 9,q) satisfies

E(0, 9, q) = 0.5, E(0,9, (?) = 0, E(+oo, 9, q) = 1

for all 9 > 0 , 0 < q < q° .

Assumption 3: For all 9 > 0, ^E _1(Z, 9,q) —> 0 (q—> 0 ) uniformly on compact subsets of the interval Z Є (0,1), and E-1(Z, 9,q) —> 9 (q —> 0 ) pointwise for all Z Є (0,1) and 9 > 0 .

Assumption 4: For all 9 > 0, the length of the interval [yi(q), y2{q)], where yi(q) inf {у Є R : E(y, 9, q) = 0} and y^{q) sup{y Є R : E(y, 9, q) — 1} , tends to 0 as q -* 0 .

3. Obtaining a system of ordinary differential equations

A method to study System (1) is well-known in the literature, and it is usually called "the linear chain trick"(see e.g. [2]). However, a direct application of this "trick"in its standard form is not suitable for our purposes, because we want Z\, as in (1), to depend on a single variable, i.e. on y\. Modifying the linear chain trick we can remove this drawback of the method.

For the sake of notational simplicity we will replace System (1) with the following scalar differential equation depending on a single response function:

x(t) = F{Z) - G(Z)x(t),

Z = Z(y),

y(t) = (9fcr)(<) {t ^ 0),

where we assume that y — yi, q = qi, 9 — 9\ Z\ — Z = E(y, 9, q).

Let us put ^

Wu{t) = I Ky(t - s)x(s)ds,

J — OC

where t ^ 0.

It is easy to see that

w\ -- —aw\ + ax and wu — aw^-i — awu, u > 2.

In what follows, we will use the following variables:

p p—v+l

y = C0X + Y, and V''—'Yh Ci+U~= 2’

3=1

(4)

(5)

(6)

i/=i

Then

x(t) = F(Z) - G{Z)x(t), v(t) = Av(t) + П(а;(і)), t> 0

(7)

where

A =

Z = E (y,9,q), У

—a a 0 ... 0

0 —a a 0

0 0 —a ... 0

0 0 ... 0 —a

, v =

/ «і \

V2

V Vp j

(8)

and

with

II(a:) := axn + Cof(Z, x)

/ Co + Cl \ / F(Z) - G'tZJa: \

C2 , f(Z,s):= 0

7T =

V Cp V 0 j

(9)

(10)

4. Stationary points

We are studying the delay system (1), which is now replaced by the equivalent system of ordinary differential equation (7).

Definition 1

A point P° is called a stationary point for System (7) with Z, = S(y*, 0j,O) (i = 1,n) if there exist a number e > 0 and points Pq , q = (<71,qn), qi G (0, e) (i = 1, ...,n) such that

• Pq is a stationary point for System (7) with Z; — T,(yi,0i,qi) (* = 1,...,n);

• Pq —> P° as q —> 0 .

It is evident that if the limit point P° does not belong to the discontinuity set of System(7) with Zi = S(j/j,0j,O), i.e. if Xi ^ 0i (i = 1, ...,n), then P° is just a usual stationary point for this system. Clearly, neither the delay operator 5R, nor the logoids Zt — S(j/,,0i,O) (i = 1 q > 0), satisfying Assumptions 1-4, influence the position of the stationary point.

Thus obtained P° is called regular stationary point [1], [4].

The case where some of the coordinates coincide with the respective thresholds is more involved. Below we provide a sufficient condition for P° to be singular stationary point [1], [4|. Proposition 1

Let Bi (i = 2, ...n) be a finite sequence consisting of 0 or 1. Assume then that

J ~ i(Zi, B2,..., Bn) - JLg^Z^B2,Bn)01

is non-zero (the derivative does not depend on Z\) and that the system

Fi(Z\, B2,..., Bn) - Gi(Zi,B2, ..., Bn)6\ = 0 Fi(Zu B2,..., Bn) - G^, B2)Bn)xi = 0 (i > 2)

(11)

(12)

with the constraints

0 < Zi < 1,

£(®i,0i,O) = Bi (i£ 2)

(13)

has a solution Z\ , x® (i ^ 2).

Then there exists a unique p-vector v° such that the point P° = (xj1,..., 2^, v°) is SSP for System (7) with Zt = 'E(yl,9l,Q) (i = 1 ,...,n). This point is independent of the choice of the operator 3? and the logoids Zj = E(j/i, %) (qi > 0, i = 1,n), satisfying Assumptions 1-4.

In a similar way, we define the notion of a stable stationary point (see e.g. [3]).

Definition 2

A stationary point P° = (x^, ...,x^,w°) for System (7) with Zt — S(yj,6>i,0) (i = l,...,n) is

called asymptotically stable if for any set of approximating stationary points Pq -> P° (q -> 0)

for System (7) with Zt = S(y*, 6hqi) (qi > 0, i = 1,n), there exist a number e > 0 such that

Pq are asymptotically stable for qi G (0, e) (i = 1, ...,n).

5. Stability results

In the non-delay case any regular stationary point is always asymptotically stable as soon as it exists. This is due to the assumptions Gi > 0, while the condition J < 0 (see (11)) gives asymptotic stability of singular stationary points (see e.g. [5]).

Including delays leads to more complicated stability conditions. We start with the case p = 1. Proposition 2

Let p = 1 and let the equation

F(Z) - G(Z)0 = 0 (14)

have a solution Z° satisfying 0 < Z° < 1.

Then the point P°(x°,y°), where x° = y° = 9, will be asymptotically stable if J < 0, and unstable if J > 0, where

J = F'{Z) - G'(Z)9 (15)

is independent of Z (as both F and G are affine).

Proposition 3

Let p — 2 and let the equation (14) have a solution Z° satisfying 0 < Z° < 1.

A. Assume that co > 0 in (2). Then the point P°(x°, y°, v°), where x° — y° = 9, vq = c29, will be asymptotically stable if J < 0, and unstable if J > 0. Moreover, assuming J < 0 there exists e > 0 such that

1. if c2 < 4coC2 , then the stationary points Pq are stable spiral points for all 0 <q<£\

2. if c2 > 4coC2 , then the stationary points Pq are stable nodes for all 0 < q < e .

B. Assume that Co = 0 in (2). Then the point P0(x°,y°,v0), where x° = y° — 9, vq = c29, has the following properties

1. If J > 0, then P° is unstable.

2. If J < 0 , ci = 0 , then P° is unstable.

3. If J < 0, ci ^ 0 and G(Z°) < aci_1(l — 2ci), then P° is unstable.

4. If J < 0, ci ^ 0 and G(Z°) > aci_1(l - 2ci), then P° is asymptotically stable (in fact,

the stationary points Pq are stable spiral points for small q > 0).

Here J is again given by (15).

The analytical formulas for n = 3,4,5 can be obtained with the help of CAS based on Mathematica.

Proposition 4

For p — 3 the Jacoby matrix is

J =

where g(q), D(q), d(q) are given by

/ -9{q) D{q)d(q) 0 0 \

a(°c +1c) —°cg(q) -a +°cD(q)d(q) a 0

a2c 0 —a a

V a3c 0 0 —a /

(16)

g(q) = G(Z(q))

D(q) = F'(Z(q)) - G\Z{q))x{q\ (17)

d(?) = H(y(?M>g)-

If

A. °c > 0, D< 0 and 9°c2 +xc(2 Jc +2c) +°c(9 + 32c - 1) > 0 or

B. °c = 0, D < 0 and a(xc —2c) +xc G(°Z) > 0, then the point P is asymptotically stable.

Вестник ТГУ, т.15, вып.2,2010_____________

Proposition 5

For р = 4 the Jacoby matrix is

J =

/ -G(°Z) D(q)d(q) 0 0 0 \

a(°c +1c) —°cg(q) -a +°cD(q)d(q) a 0 0

a2c 0 —a a 0

a3c 0 0 —a a

V OL ^c 0 0 0 —a /

(18)

where g(q), D(q), d(q) are given by (17).

If

A. °c > 0, D < 0 , 20 °c2 +lc (3 lc +2c) +°c (15 lc + 2 2c —3c) > 0 and 80 °c3 + 8 °c? (15 h + 6 2c +

2 3c - 2) +xc (91c2 +2c (2 2c +3c) +1c (9 2c + 3 3c - 1)) +°c (57 *c2 + 4 2c2 -3c2 + 4 (10 2c + 3 3c - 2)) >0

or

B. °c = 0, D < 0 , a{lc -2c) +1cG(°Z) > 0 and

9 +2c (22c +3c) +xc(92c + 33c — 1) > 0,

then the point P is asymptotically stable.

Proposition 6

For p = 5 the Jacoby matrix is

V

—G(°Z) a(°c+1c) —°cg(q) D(q)d(q) 0 0 0 0

-a +°cD(q)d(q) a 0 0 0

a 2c 0 —a a 0 0

a 3c 0 0 —a a 0

a4c 0 0 0 —a a

a 5c 0 0 0 0 —a

(19)

where g(q), D(q), d(q) are given by (17).

If

A. °c > 0, D < 0 , 40 °c2 +xc (4 +2c) +°c (24 lc + 2 2c —3c) > 0,

275°c3 +°c (139 V + (2 2c -3c)(3 2c +3c) +xc (642c - 2% - 10^ + 1)) +xc (20 :c2 + %(32c +3c) + fc (15 2c + 2 3c —4c)) + 5 °c2 (66 lc + 13 2c — 4 3c — 5 4c + 1) >0 and

1375 °c4 + 50 °c2 (55 Jc+23 2c+9?; + 3 4c - 7) +°c2 (2015 h2 + 225 2c2 + 30 - 45 3c2 -1 + 5 2c (13 3c -24c — 6) + 10^ - 703c4c - 25 V + 10 Jc(1572c + 573c + 184c - 35)) +2c(80 rc3 + 8xc2 (152c + 63c + 2^—2) +2c (9 2c2 +3c (2 3c +4c) +2c (9 3c + 3 — 1)) +1c (57 2c2 + 4 3c2 + 4 % (10 3c + 3 — 2))) +

°c(656+ 2\? (3742c + 1403c + 47*e - 64) + (22c -3c)(92c? +3c(23c +4c) +2c(93c + 3*c - 1)) + (2312c2 + 2c (123 3c + 344c - 36) + 2 (43c - 43c2 +4c - 113c^ - 5 \?))) > 0 or

B. °c = 0, D < 0, a(}c —2c) +1cG(°Z) > 0, 20 V2 +2c (32c +3c) +xc (152c + 23c —4c) > 0 and 80 h3 + 8 V2 (15 2c + 6 3c + 2 ^ - 2) +2c(9 2c2 +3c (2 3c +^c) +2c (9 b + 3 *c - I)) +Jc (57 2c2 + 4 3c2 -

iP 4- 4 2c (10 3c + 3 y — 2)) > 0, G(°Z) > 0, then the point P is asymptotically stable.

Finally, we remark that a more sophisticated asymptotic analysis can be performed on the basis of parallel computer algebra routines.

REFERENCES

[1] L. Glass and S. A. Kauffman . The logical analysis of continuous, non-linear biochemical control networks, J. Theor. Biol., v. 39 (1973), 103-129.

[2] N. McDonald. Time lags in biological models, Lect. Notes in Biomathematics, Springer-Verlag, Berlin-Heidelberg-New York, 1978, 217 p.

[3] T. Mestl, E. Plahte, and S. W. Omholt, A mathematical framework for describing and analysing gene regulatory networks, J. Theor. Biol., v. 176 (1995), 291-300.

[4] E. Plahte, T. Mestl, and S. W. Omholt, Global analysis of steady points for systems of

differential equations with sigmoid interactions, Dynam. Stabil. Syst., v. 9, no. 4 (1994), 275-291.

Вестник ТГУ, т.15, вып.2, 2010

[5] Е. Plahte, Т. Mestl, and S. W. Omholt, A methodological basis for description and analysis of systems with complex switch-like interactions, J. Math. Biol., v. 36 (1998), 321-348.

[6] A. PONOSOV, Gene regulatory networks and delay differential equations. Special issue of Electronic J. Diff. Eq., v. 12 (2004), pp. 117-141.

Acknowledgment

The authors would like to thank The Norwegian Programme for Development, Research and Education (NUFU) for partial financial support.

It is received on March, 21, 2010

Шлыкова И., Поносов А., Непомнящих Ю., Шиндяпин А. Анализ устойчивости в окрестностях сингулярных стационарных точек и компьютерные алгебраические системы. Предлагается метод изучения асимптотических свойств решений систем дифференциальных уравнений с распределённым запаздыванием и булевыми нелинейностями (функции скачков). Подобные системы возникают во многих приложениях, предметом этой работы являются системы, описывающие регулируемые генные сети. Спецификой является анализ устойчивости стационарных точек, принадлежащих множествам разрыва систем (пороговым состояниям). Анализируется устойчивость таких точек в случае запаздывания по временной переменной. Основным методом является замена ступенчатых функций их гладкими аппроксимациями, так называемыми «логоидами», и изучение получаемых таким образом гладких систем.

Ключевые слова: регулируемые генные сети, ступенчатые функции, устойчивость разрывных систем.