Models of genetic networks with given properties

O. KOZLOVSKA1, F. SADYRBAEV1,2

1Department of Natural Sciences and Mathematics

Daugavpils University

Parades street 1, Daugavpils

LATVIA,

2Institute of Mathematics and Computer Science,

University of Latvia

Rainis boulevard 29, Riga

LATVIA

Abstract: A multi-parameter system of ordinary diﬀerential equations, modelling genetic networks, is

considered. Attractors of this system correspond to future states of a network. Suﬃcient conditions for

the non-existence of stable critical points are given. Due to the special structure of the system, attractors

must exist. Therefore the existence of more complicated attractors was expected. Several examples are

considered, conﬁrming this conclusion.

Key-Words: Mathematical model, dynamical system, attractors, genetic regulatory networks

Received: March 17, 2021. Revised: February 18, 2022. Accepted: March 19, 2022. Published: April 20, 2022.

1 Introduction

We consider systems of ordinary diﬀerential equa-

tions of the form











dt =1

1+e−µ1(w11 x+w12 y+w13 z−θ1)−v1x,

dt =1

1+e−µ2(w21 x+w22 y+w23 z−θ2)−v2y,

dt =1

1+e−µ3(w31 x+w32 y+w33 z−θ3)−v3z

(1)

Such systems arise in the theory of complex net-

works, such as telecommunication ones [4], ge-

netic networks [3], [2], [7], [1], [6], [11], [12], neu-

ronal networks [9], [5]. More information can be

found in the review papers [3], [8] and others. In

models of genetic regulatory networks (GRN in

short) nonlinear functions on the right side de-

scribe interrelation between genes. They commu-

nicate by protein expression, sending messages to

each other. In such a way a collective reaction of

a network is formed. In the system (1) the func-

tion f(z) = 1

1+e−µz is sigmoidal. Sigmoidal func-

tions are monotonically increasing from zero to

unity. These functions are most suitable for mod-

eling purposes. The parameters µiare positive.

The linear part in (1) describes the natural decay

of a network if there are no connections between

genes. The parameters θcan be any numbers and

in many cases they are adjustable. The coeﬃ-

cients wij are the entries of the regulatory matrix

Wthat describes links between genes. A positive

element wij means activation the i-th gene by j-

th one. Negative coeﬃcient means repression (in-

hibition), zero means no relation. The change of

one entry of matrix Wleads to rearrangement

of the whole network. The collective response of

GRN to changes in environment can be developed

quickly. The current state of GRN is represented

by the state vector X(t) = (x1, . . . , xn(t)).The

dimension nin system (1) is three. Mathemat-

ically X(t) is the solution vector of system (1).

Properties of system (1) are such that it has an

invariant set in the phase space R3.This set is a

parallelepiped Q3={0< xi<1/vi, i = 1,2,3}.

Due to the properties of sigmoidal functions the

vector ﬁeld deﬁned by system (1), is directed in-

wards of the region Q3.The trajectories X(t) are

therefore trapped in Q3.The attractors of system

(1) must exist in Q3.Future states of a network

are heavily dependent on the number, location

and properties of attractors. The structure and

properties of attractors for system (1) are to be

studied.

Equilibria of system (1) always exist and are

located in the domain Q3.Equilibria are solutions

WSEAS TRANSACTIONS on COMPUTER RESEARCH

DOI: 10.37394/232018.2022.10.6

O. Kozlovska, F. Sadyrbaev

E-ISSN: 2415-1521

Volume 10, 2022

of the system











v1x=1

1+e−µ1(w11 x+w12 y+w13 z−θ1),

v2y=1

1+e−µ2(w21 x+w22 y+w23 z−θ2),

v3z=1

1+e−µ3(w31 x+w32 y+w33 z−θ3).

(2)

The system (2) consist of three equations which

deﬁne the nullclines. Three nullclines can inter-

sect only in Q3.Indeed, suppose that the critical

point (x∗, y∗, z∗) locates outside Q3.Then some

coordinate, say x∗,either is negative or greater

than 1/v1.This contradicts the properties of the

sigmoidal function in the right side (recall that

the range of values of a sigmoidal function is

(0,1)). The number of intersections is ﬁnite, but

cannot be zero. A unique cross-point is possible,

as examples show.

Our intent in this paper is to clarify the follow-

ing situation. Imagine that there is a unique crit-

ical point (an equilibrium), which is not attrac-

tive. The suﬃcient condition for non-attractivity

is positivity of a real part of at least one char-

acteristic number λ. The characteristic numbers

can be computed for any equilibrium using the

standard analysis of the linearized equation. We

will provide the necessary information in the next

section. If a unique critical point is not attractive,

there are other attractors.

We pass to description of our results in this

direction. We provide the suﬃcient conditions,

in terms of the coeﬃcients of Wand parame-

ters of the system (1), for a critical point to be

non-attractive. We construct then the examples,

where our conditions are applicable. In these ex-

amples the critical point is unique. Therefore an

attractor of more complicated structure should

exist. We have discovered them. To be more spe-

ciﬁc, our intent is to provide conditions for ex-

istence of equilibria with the characteristic num-

bers λ1∈R, λ2,3=α±βi, where α > 0, i =√−1.

2 Problem Formulation

Consider system (1), where v1=v2=v3= 1,for

simplicity.

The critical points are solutions of the system











x=1

1+e−µ1(w11 x+w12 y+w13 z−θ1)

y=1

1+e−µ2(w21 x+w22 y+w23 z−θ2),

z=1

1+e−µ3(w31 x+w32 y+w33 z−θ3).

(3)

The linearized system at a critical point

(x∗, y∗, z∗) is











1=−u1+µ1w11g1u1

+µ1w12g1u2+µ1w13g1u3,

2=−u2+µ2w21g2u1

+µ2w22g2u2+µ2w23g2u3,

3=−u3+µ3w31g3u1

+µ3w32g3u2+µ3w33g3u3,

(4)

where

g1=e−µ1(w11 x∗+w12 y∗+w13 z∗−θ1)

[1 +e−µ1(w11 x∗+w12 y∗+w13 z∗−θ1)]2,

g2=e−µ2(w21 x∗+w22 y∗+w23 z∗−θ2)

[1 +e−µ2(w21 x∗+w22 y∗+w23 z∗−θ2)]2,

g3=e−µ3(wn1x∗+wn2y∗+w33 z∗−θ3)

[1 +e−µ3(w31 x∗+w32 y∗+w33 z∗−θ3)]2.

Properties of a critical point (x∗, y∗, z∗) are de-

scribed by the three numbers (they are called the

characteristic numbers) λ1, λ2, λ3,which can be

found from the chracteristic equation. One has

det(A−λI) = 0,(5)

where Ais 3 ×3 matrix of the coeﬃcients of the

system (4), Iis the unity 3 ×3 matrix.

The equation (5) is complicated, since it re-

quires computation of the partial derivatives gi

at the critical point (x∗, y∗, z∗).It can be sim-

pliﬁed, however, in this way. Observe, from (3),

that

e−µ1(w11 x∗+w12 y∗+w13 z∗−θ1)=1

x∗−1,

e−µ2(w21 x∗+w22 y∗+w23 z∗−θ2)=1

y∗−1,

e−µ3(w31 x∗+w32 y∗+w33 z∗−θ3)=1

z∗−1.

(6)

Then g1=

x∗−1

x∗2

=(1 −x∗)x∗and, similarly,

g2= (1 −y∗)y∗, g3= (1 −z∗)z∗.

The entries aij of the matrix A−λI, needed for

construction of the characteristic equation, are:

a11 =µ1w11(1 −x∗)x∗−1−λ,

a12 =µ1w12(1 −x∗)x∗,

WSEAS TRANSACTIONS on COMPUTER RESEARCH

DOI: 10.37394/232018.2022.10.6

O. Kozlovska, F. Sadyrbaev

E-ISSN: 2415-1521

Volume 10, 2022

a13 =µ1w13(1 −x∗)x∗,

a21 =µ2w21(1 −y∗)y∗,

a22 =µ2w22(1 −y∗)y∗−1−λ,

a23 =µ2w23(1 −y∗)y∗,

a31 =µ3w31(1 −z∗)z∗,

a32 =µ3w32(1 −z∗)z∗,

a33 =µ3w33(1 −z∗)z∗−1−λ.

Choice of θ-s.The coordinates of a critical

point (x∗, y∗, z∗) are solutions of the system (3).

Observe, that if x∗=y∗=z∗= 0.5,then for θi

such that

w11 +w12 +w13 = 2 θ1,

w21 +w22 +w23 = 2 θ2,

w31 +w32 +w33 = 2 θ3,

(7)

the numbers x∗=y∗=z∗= 0.5 solve the system

(3).

Assumption θ.We assume that θisatisfy

(7).

Conclusion Θ.The entries aij of the matrix

A−λI, needed for construction of the character-

istic equation, become:

a11 =µ1w110.25 −1−λ,

a12 =µ1w120.25,

a13 =µ1w130.25,

a21 =µ2w210.25,

a22 =µ2w220.25 −1−λ,

a23 =µ2w230.25,

a31 =µ3w310.25,

a32 =µ3w320.25,

a33 =µ3w330.25 −1−λ.

Choice of µ-s. For easier calculations, we

would like to make the products µi0.25 in the

above elements equal to unity. Therefore we

choose µi= 4 for i= 1,2,3.

Assumption µ.We assume that µi= 4 for

i= 1,2,3.

Conclusion µ.The entries aij of the matrix

A−λI, needed for construction of the character-

istic equation, become:

a11 =w11 −1−λ,

a12 =w12,

a13 =w13,

a21 =w21,

a22 =w22 −1−λ,

a23 =w23,

a31 =w31,

a32 =w32,

a33 =w33 −1−λ.

2.1 Characteristic equation

We are now in a position to construct the

characteristic equation for the critical point

(0.5,0.5,0.5) assuming that Assumption θand

Assumption µhold. Computing the determinant

of the matrix A−λI, we obtain the equation

Λ3−(w11 +w22 +w33)Λ2

−(w21w12 −w11w22 +w31w13

+w32w23 −w11w33 −w22w33)Λ

−(−w31w22w13 +w21w32w13

+w31w12w23 −w11w32w23

−w21w12w33 +w11w22w33) = 0,

(8)

Λ = λ+ 1.(9)

Rewrite (8) in a compact form

Λ3+PΛ2+QΛ + R= 0,(10)

where

P=−(w11 +w22 +w33),(11)

Q=−(w21w12 −w11w22 +w31w13

+w32w23 −w11w33 −w22w33),(12)

R=−(−w31w22w13 +w21w32w13 +w31w12w23

−w11w32w23 −w21w12w33 +w11w22w33).

(13)

Equation (10) for the variable λis

λ3+(3+P)λ2+(3+2P+Q)λ+(1+P+Q+R) = 0.

(14)

2.2 Problem

We can now formulate the problem. Up to now

we have made the choice of the parameters vi, θi

and µi.The critical point under consideration is

placed to the location (0.5,0.5,0.5) by the choice

of θi.

Our plan is the following:

1) we make the point (0.5,0.5,0.5) non-

attractive by the appropriate choice of the

entries in the matrix W;

2) under the assumption that this point is

unique we are led to the conclusion, that there

are no attractive equilibria in the cube Q3;

3) Then there must exist other attractors.

3 Problem Solution

We assume that there is a single critical point and

it is placed in the middle point by an appropri-

ate choice of the parameters θ. We wish to make

this point non-attractive by the choice of other

parameters and matrix W.

WSEAS TRANSACTIONS on COMPUTER RESEARCH

DOI: 10.37394/232018.2022.10.6

O. Kozlovska, F. Sadyrbaev

E-ISSN: 2415-1521

Volume 10, 2022

3.1 Routh - Gurwitz criterion and

related material

Let us recall the Routh-Gurwitz conditions for

the third order equations. Consider

a3s3+a2s2+a1s+a0= 0.(15)

Our trusted source of information is [10].

Consider the Table:

s3a3a1

s2a2a0

sa2a1−a0a3

s0a0

(16)

This table is called Routh array.

Proposition 3.1. The number of roots in the

open right half-plane is equal to the number of

sign changes in the ﬁrst column of Routh array.

Proposition 3.2. If 0 appears in the ﬁrst col-

umn in nonzero row in Routh array, replace it

with a small positive number. In this case, poly-

nomial has some roots in the open right half-plane

(RHP).

Proposition 3.3. If zero row appears in Routh

array, polynomial has roots either on the imagi-

nary axis or in RHP.

Proposition 3.4. If all zeros are in the left half

plane (LHP), then all akare of the same sign.

Corollary 1. If there is sign change in akthen

some zero is not in the LHP.

Theorem 1. If (3 + P)or (3 + 2P+Q)or (1 +

P+Q+R)is negative in (14), then some λis

not in the LHP.

Example 1. 3 + P < 0.

Take µi= 4,3−(w11 +w22 +w33)<0.All

others wij are taken randomly. To be speciﬁc, let

w11 = 1, w22 = 2, w33 = 3.Let wij = +1 for i+j

even, wij =−1 for i+jodd. So the matrix

W=Ã1−1 1

−1 2 −1

1−1 3 !.(17)

is suitable. Choose θias recommended, so θ1=

0.5, θ2= 0.0, θ3= 1.5.The characteristic equa-

tions is

λ3+ (3 + P)λ2+(3+2P+Q)λ+ (1 + P+Q

+R) = λ3−3λ2−λ+ 1 = 0.

(18)

The Routh array is

s31−1

s2−3 1

s01

(19)

There are two sign changes in the ﬁrst column of

the Routh array.

Solving (18) with respect to λ, obtain

λ1=−0.675131, λ2= 0.460811, λ3= 3.21432.

The type of the equilibrium at (0.5,0.5,0.5) is a

3D saddle point. There are other critical points,

so the existence of attractors other than the sta-

ble equilibria is not observed.

Example 2. 3+2P+Q < 0.

Let the matrix be

W=Ã3 5 2

−1 1 0

0 1 −2!.(20)

The parameters µiand θiare chosen as usu-

ally. After simple calculations obtain 3 + P= 1,

3 + 2P+Q=−1,1 + P+Q+R= 17.The

characteristic numbers are λ1=−3.1, λ2,3=

1.05 ±2.1i.The critical point at (0.5,0.5,0.5) is

non-attractive.

Example 3. 1 + P+Q+R < 0.

Choose six entries in the regulatory matrix and

let the remaining three be arbitrary. Suppose

wii = 0 for i= 1,2,3.Let w13 =−1, w21 =−1,

w32 =−1.So the matrix Wis of the form

W=Ã0w12 −1

−1 0 w23

w31 −1 0 !.(21)

Then

P= 0, Q =w12 +w31 +w23, R = 1−w31w12w23,

1+P+Q+R= 1+w12+w31 +w23+1−w31w12w23

and the condition 1 + P+Q+R < 0 takes the

form

w12 +w31 +w23 < w31w12w23 −2.(22)

The region is visualized in two pictures Fig. 1

and Fig. 2.

The elements w12 =−2, w23 =−2, w31 = 1

are admissible.

For the matrix

W=Ã0−2−1

−1 0 −2

1−1 0 !,(23)

WSEAS TRANSACTIONS on COMPUTER RESEARCH

DOI: 10.37394/232018.2022.10.6

O. Kozlovska, F. Sadyrbaev

E-ISSN: 2415-1521

Volume 10, 2022

Figure 1: First view.

Figure 2: Second view.

parameter µi=4, θi= (wi1+wi2+wi3)/2

the characteristic equation for the central point

(0.5,0.5,0.5) is

λ3+ 3λ2−0λ−5 = 0.(24)

The Routh array is

s31 0

s23−5

s0−5

(25)

There is one sign change in the ﬁrst column.

The characteristic numbers for the point

(0.5,0.5,0.5) are

λ1= 1.1038, λ2,3=−2.0519 ±0.565236i.

The type of this critical point is an attractive

focus and repelling in the remaining dimension.

The point is not attractive, but there are other

critical points and the existence of an attractor

other than a stable equilibrium, cannot be guar-

anteed.

In the above examples we succeded in assign-

ing the critical point (0.5,0.5,0.5) the status of

non-attracting point. There are, however, other

critical points and some of them (or all of them)

are attractive. In the next section we provide ex-

amples, where assumtions 1 and 2 of subsection

2.2 together are satisﬁed.

4 More examples

I. Consider system (1), where µi= 4, vi= 1,and

the regulatory matrix

W=Ã1−a−1

−1 1 −a

−a−1 1 !.(26)

Suppose a≥0.If a= 0,the matrix becomes

W=Ã1 0 −1

−1 1 0

0−1 1 !.(27)

and it contains inhibitory cycle (three −1) and

self-activation (three 1). The corresponding sys-

tem (1) is known to have periodic solutions.

Consider system (1) with the matrix (28).

Suppose that the parameters θiare chosen ap-

proprietly.

Computational experiments show that for a=

0.1,0.2,0.3,0.4, ..., 0.8 there exists a periodic so-

lution. There exists a single critical point at

(0.5,0.5,0.5).

For a= 0.8 the values (3 + P),(3 + 2P+

Q),(1 + P+Q+R),mentioned in Theorem 1,

respectively are 0,−2.4,1.512.Therefore a single

critical point has the real part in RHP.

The periodic solution exists also, depicted in

the Figure 3.

0.0

0.5

1.0

0.0

0.5

1.0

0.0

0.5

1.0

Figure 3: Stable periodic solution for a= 0.8

II. Consider Example 2 in the previous sec-

tion. We have shown that the critical point at

the center is not attractive. Visual inspection of

WSEAS TRANSACTIONS on COMPUTER RESEARCH

DOI: 10.37394/232018.2022.10.6

O. Kozlovska, F. Sadyrbaev

E-ISSN: 2415-1521

Volume 10, 2022

nullclines conﬁrm that this point is unique. It has

the following characteristic numbers λ1=−3.10,

λ2,3= 1.05 ±2.10i.

Figure 4: The nullclines in Example II.

Therefore another attractor is expected. It ex-

ists in the form of stable periodic solution.

0.0

0.5

1.0

0.0

0.5

1.0

0.0

0.5

1.0

Figure 5: Stable periodic solution in Example

II.

III. Consider system (1) with the previous

choice of parameters and the regulatory matrix

W=Ã3 5 2

−110

0 1 1 !.(28)

The values (3+P),(3+ 2P+Q),(1+P+Q+R),

mentioned in Theorem 1, are −2,5,2.Therefore a

single critical point has the real part in RHP. The

craracteristic numbers are λ1=−0.34, λ2,3=

1.17 ±2.11i.An attractor exists.

Figure 6: Stable periodic solution in Example

III.

5 Conclusion

Using mathematical models can make the process

of studying GRN easier and more eﬀective. Given

experimental data of some particular GRN, the

mathematical model should be adjusted by the

appropriate change of parameters. Afterward,

the model can be studied using standard math-

ematical analysis tools. An inverse problem is

to take an existed model and try to predict the

properties of a real network, studying its model.

In our research, we did something in this spirit.

We try to select the characteristics of a network

(respectively, the elements of the regulatory ma-

trix W) in such a way, that the existence of sta-

ble equilibria is not allowed. This can reveal

more complicated attractive sets on the phase

space. We have obtained such conditions and

tested them. In our examples, which were con-

structed, following the recommendations of The-

orem 1, no stable equilibria can appear. Instead,

we got stable periodic attractors. In perspective,

considering models with several critical points of

non-attractive nature, we hope to get attractors

of the next level of complicacy. Recall, that the

popular Lorenz attractor is based on three equi-

libria, two of them are non-attractive saddle fo-

cuses and one is a saddle point.

References:

[1] E. Brokan, F. Sadyrbaev. Remarks on GRN

type systems. 4open, vol. 3 (2020), arti-

cle number 8. https://doi.org/10.1051/

fopen/2020009

WSEAS TRANSACTIONS on COMPUTER RESEARCH

DOI: 10.37394/232018.2022.10.6

O. Kozlovska, F. Sadyrbaev

E-ISSN: 2415-1521

Volume 10, 2022

[2] C. Furusawa, K. Kaneko. A generic mech-

anism for adaptive growth rate regula-

tion. PLoS Comput Biol 4(1), 2008: e3.

doi:10.1371/journal.pcbi.0040003

[3] H.D. Jong. Modeling and Simulation of Ge-

netic Regulatory Systems: A Literature Re-

view, J. Comput Biol. 2002;9(1):67-103, DOI:

10.1089/10665270252833208

[4] Y. Koizumi et al. Adaptive Virtual Network

Topology Control Based on Attractor Selec-

tion. J. of Lightwave Technology, (ISSN :0733-

8724), Vol.28 (06/2010), Issue 11, pp. 1720-

1731 DOI:10.1109/JLT.2010.2048412

[5] V.W. Noonburg Diﬀerential Equations: From

Calculus to Dynamical Systems, Providence,

Rhode Island: MAA Press, 2019, 2nd edition.

[6] F. Sadyrbaev, S. Atslega, E. Brokan. (2020)

Dynamical Models of Interrelation in a Class

of Artiﬁcial Networks. In: Pinelas S., Graef

J.R., Hilger S., Kloeden P., Schinas C. (eds)

Diﬀerential and Diﬀerence Equations with

Applications. ICDDEA 2019. Springer Pro-

ceedings in Mathematics & Statistics, vol

333. Springer, Cham. https://doi.org/10.

1007/978-3-030-56323-3_18

[7] F. Sadyrbaev, V. Sengileyev, A. Sil-

vans. On Coexistence of Inhibition and

Activation in Genetic Regulatory Net-

works. 19 Intern.Confer.Numer.Analys.

and Appl.Mathematics, Rhodes,

Greece, 20-26 September 2021, To ap-

pear in AIP Conference Proceedings.

https://aip.scitation.org/journal/apc

[8] N. Vijesh, S.K. Chakrabarti, J. Sreekumar.

Modeling of gene regulatory networks: A re-

view, J. Biomedical Science and Engineering,

6:223-231, 2013.

[9] H.R. Wilson, J.D. Cowan. Excitatory and in-

hibitory interactions in localized populations

of model neurons. Biophys J., vol 12 (1), 1972,

pp. 1-24.

[10] https://www.egr.msu.edu/classes/

me451/jchoi/2012/notes/ME451_L10_

RouthHurwitz.pdf

[11] F. Sadyrbaev, I. Samuilik, V. Sengileyev. On

Modelling of Genetic Regulatory Networks.

WSEAS Transactions on Electronics, vol. 12,

pp. 73-80, 2021

[12] I. Samuilik, F. Sadyrbaev. Mathematical

Modelling of Leukemia Treatment. WSEAS

Transactions on Computers, vol. 20, pp. 274-

281, 2021.

Contribution of individual

authors to the creation of a

scientiﬁc article (ghostwriting

policy)

Author Contributions: Please, indicate

the role and the contribution of each

author:

Example

Both authors have contributed equally to cre-

ation of this article.

Follow: www.wseas.org/multimedia/contributor-

role-instruction.pdf

Sources of funding for research

presented in a scientiﬁc article

or scientiﬁc article itself

Report potential sources of funding if there

is any

Creative Commons Attribution

License 4.0 (Attribution 4.0

International, CC BY 4.0)

This article is published under the terms of the

Creative Commons Attribution License 4.0

https://creativecommons.org/licenses/by/4.0/deed.en US

WSEAS TRANSACTIONS on COMPUTER RESEARCH

DOI: 10.37394/232018.2022.10.6

O. Kozlovska, F. Sadyrbaev

E-ISSN: 2415-1521

Volume 10, 2022