INNA SAMUILIK
Department of Natural Sciences and Mathematics
Daugavpils University
Parades street1
LATVIA
Abstract: Mathematical modeling is a universal tool for the study of complex systems.
In this paper formulas for characteristic numbers of critical points for the systems of
order four (4D) are considered. We show how an unstable focus-focus can appear in a
four-dimensional system. Projections of 4D trajectories on two-dimensional and three-
dimensional subspaces are shown. In the considered four-dimensional system the logistic
function is used. The research aims to investigate the four-dimensional system, find
critical points of the system, calculate the characteristic numbers, and calculate Lyapunov
exponents.
Key-Words: mathematical modeling, nullclines, Lyapunov exponents, periodic solutions,
linearized system, critical points
Received: April 29, 2021. Revised: October 13, 2022. Accepted: November 15, 2022. Published: December 20, 2022.
1 Introduction
Every natural science consists of three
parts: empirical, theoretical, and math-
ematical. The empirical part contains
factual information obtained in experi-
ments and observations. The theoretical
part develops theoretical concepts. The
mathematical part constructs mathema-
tical models that serve to test the basic
theoretical concepts provide methods for
the primary processing of experimental
data so that they can be compared with
the results of the models, and develops
methods for planning an experiment in
such a way that, with a small expenditure
of effort, it is possible to obtain sufficiently
reliable data from experiments, [1]. Mathe-
matical models are routinely used in the
physical and engineering sciences to help
understand complex systems and optimize
industrial processes. There are numerous
examples of the fruitful application of
mathematical principles to problems in cell
and molecular biology, and recent years
have seen increasing interest in applying
quantitative techniques to problems in
biotechnology, [3]. The main problem in
the mathematical modeling of a dynamic
system is to develop a model and then to
determine dependencies and coefficients in
the equations used in developing the model,
[7]. To quote the statistician Dr. George
E. P. Box (1919-2013): “Essentially, all
models are wrong, but some are useful, [8].”
Gene regulatory networks are structurally
represented by spatially located objects,
consisting of occurring and hundreds of
elements of different natures and com-
plexity. The most important property
of the gene regulatory networks is the
ability to change the state in response to
changes in the conditions of the external
and internal environment, [2]. A variety of
approaches have been used to model the
gene regulatory networks: Graph method,
differential, statistical equations, and more.
Mathematical Modeling of Four-dimensional Genetic Regulatory
Networks Using a Logistic Function
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.15
Inna Samuilik
E-ISSN: 2415-1521
112
Volume 10, 2022
Consider the matrix
W=
1 1 00110
1100110
0110010
1111111
1111110
0100110
1000000
(1)
The graph with the regulatory matrix (1) is
considered in Figure 1.
Figure 1: The graph with the regulatory mat-
rix (1).
2 Four-dimensional sys-
tems
We consider systems of ordinary differential
equations of the form
dx1
dt =
1
1+e−µ1(w11x1+w12 x2+w13x3+w14 x4−θ1)−v1x1,
dx2
dt =
1
1+e−µ2(w21x1+w22 x2+w23x3+w24 x4−θ2)−v2x2,
dx3
dt =
1
1+e−µ3(w31x1+w32 x2+w33x3+w34 x4−θ2)−v3x3,
dx4
dt =
1
1+e−µ4(w41x1+w42 x2+w43x3+w44 x4−θ4)−v4x4.
(2)
Such systems arise in the theory of complex
networks, such as genetic networks, [4], [5],
[6],[7], [9], telecommunications networks,
[10], neuronal networks, [11], and more.
The greater the number of equations in
the system, the closer the model is to a
realistic gene regulatory network. But
even the three-dimensional system contains
many parameters, which makes the study
not trivial. In this paper, we consider
the four-dimensional system. This system
consists of 24 parameters.
The system (2) consists of four equations
that define the nullclines. The nullclines
are given by
v1x1=1
1+e−µ1(w11x1+w12 x2+w13x3+w14 x4−θ1),
v2x2=1
1+e−µ2(w21x1+w22 x2+w23x3+w24 x4−θ2),
v3x3=1
1+e−µ3(w31x1+w32 x2+w33x3+w34 x4−θ2),
v4x4=1
1+e−µ4(w41x1+w42 x2+w43x3+w44 x4−θ4).
(3)
Critical points are solutions of the system
(3).
2.1 Linearized system
The linearized system for critical point
(x∗
1, x∗
2, x∗
3, x∗
4) is
u0
1=−v1u1+µ1w11g1u1+µ1w12g1u2+
µ1w13g1u3+µ1w14g1u4,
u0
2=−v2u2+µ2w21g2u1+µ2w22g2u2+
µ2w23g2u3+µ2w24g2u4,
u0
3=−v3u3+µ3w31g3u1+µ3w32g3u2+
µ3w33g3u3+µ3w34g3u4,
u0
4=−v4u4+µ4w41g4u1+µ4w42g4u2+
µ4w34g4u3+µ4w44g4u4,
where
g1=e−µ1(w11x∗
1+w12x∗
2+w13x∗
3+w14x∗
4−θ1)
[1 +e−µ1(w11x∗
1+w12x∗
2+w13x∗
3+w14x∗
4−θ1)]2,
g2=e−µ2(w21x∗
1+w22x∗
2+w23x∗
3+w24x∗
4−θ2)
[1 +e−µ2(w21x∗
1+w22x∗
2+w23x∗
3+w24x∗
4−θ2)]2,
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.15
Inna Samuilik
E-ISSN: 2415-1521
113
Volume 10, 2022
g3=e−µ3(w31x∗
1+w32x∗
2+w33x∗
3+w34x∗
4−θ3)
[1 + e−µ3(w31x∗
1+w32x∗
2+w33x∗
3+w34x∗
4−θ3)]2,
g4=e−µ4(w41x∗
1+w42x∗
2+w43x∗
3+w44x∗
4−θ4)
[1 +e−µ4(w41x∗
1+w42x∗
2+w43x∗
3+w44x∗
4−θ4)]2.
Properties of a critical point (x∗
1, x∗
2, x∗
3, x∗
4)
are described by the four numbers (they
are called the characteristic numbers)
λ1;λ2;λ3;λ4which can be found from the
chracteristic equation, [9].
The characteristic equation is
λ4+Aλ3+Bλ2+Mλ +L= 0,(4)
where
A= (v1+v2+v3+v4)−g1w11µ1−
g2w22µ2−g3w33µ3−g4µ4w44,
B=v3v4−g1v3w11µ1−g1v4w11µ1−
g2v3w22µ2−g2v4w22µ2−g1g2w21w12µ1µ2+
g1g2w11w22µ1µ2−g3v4w33µ3−
g1g3w31w13µ1µ3+g1g3w11w33µ1µ3−
g2g3w32w23µ2µ3+g2g3w22w33µ2µ3−
g1g4w41µ1µ4w14 −g2g4w42µ2µ4w24−
g3g4w43µ3µ4w34 −g4v3µ4w44+
g1g4w11µ1µ4w44 +g2g4w22µ2µ4w44+
g3g4w33µ3µ4w44+
v2(v3+v4−g1w11µ1−g3w33µ3−g4µ4w44)+
v1(v2+v3+v4−g2w22µ2−g3w33µ3−g4µ4w44),
M=−g1v3v4w11µ1−g2v3v4w22µ2−
g1g2v3w21w12µ1µ2−g1g2v4w21w12µ1µ2+
g1g2v3w11w22µ1µ2+g1g2v4w11w22µ1µ2−
g1g3v4w31w13µ1µ3+g1g3v4w11w33µ1µ3−
g2g3v4w32w23µ2µ3+g2g3v4w22w33µ2µ3+
g1g2g3w31w22w13µ1µ2µ3−
g1g2g3w21w32w13µ1µ2µ3−
g1g2g3w31w12w23µ1µ2µ3+
g1g2g3w11w32w23µ1µ2µ3+
g1g2g3w21w12w33µ1µ2µ3−
g1g2g3w11w22w33µ1µ2µ3−
g1g4v3w41µ1µ4w14+
g1g2g4w41w22µ1µ2µ4w14−
g1g2g4w21w42µ1µ2µ4w14+
g1g3g4w41w33µ1µ3µ4w14−
g1g3g4w31w43µ1µ3µ4w14−
g2g4v3w42µ2µ4w24−
g1g2g4w41w12µ1µ2µ4w24+
g1g2g4w11w42µ1µ2µ4w24+
g2g3g4w42w33µ2µ3µ4w24−g2g3g4w32w43µ2µ3µ4w24−
g1g3g4w41w13µ1µ3µ4w34+g1g3g4w11w43µ1µ3µ4w34−
g2g3g4w42w23µ2µ3µ4w34+g2g3g4w22w43µ2µ3µ4w34+
g1g4v3w11µ1µ4w44 +g2g4v3w22µ2µ4w44+
g1g2g4w21w12µ1µ2µ4w44−g1g2g4w11w22µ1µ2µ4w44+
g1g3g4w31w13µ1µ3µ4w44−g1g3g4w11w33µ1µ3µ4w44+
g2g3g4w32w23µ2µ3µ4w44−g2g3g4w22w33µ2µ3µ4w44+
v1(v3v4−g2v3w22µ2−g2v4w22µ2−g3v4w33µ3−
g2g3w32w23µ2µ3+g2g3w22w33µ2µ3−
g2g4w42µ2µ4w24−
g3g4w43µ3µ4w34 −g4v3µ4w44+
g2g4w22µ2µ4w44+
g3g4w33µ3µ4w44+
v2(v3+v4−g3w33µ3−g4µ4w44))+
v2(v3(v4−g1w11µ1−g4µ4w44)−
g1µ1(v4w11 +g3w31w13µ3−
g3w11w33µ3+g4w41µ4w14 −g4w11µ4w44)−
g3µ3(v4w33 +g4w43µ4w34 −g4w33µ4w44)),
L=v1(v2(v3(v4−g4µ4w44)−
g3µ3(v4w33 +g4w43µ4w34 −g4w33µ4w44))−
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.15
Inna Samuilik
E-ISSN: 2415-1521
114
Volume 10, 2022
g2µ2(v3(v4w22 +g4µ4(w42w24 −w22w44))+
g3µ3(v4(w32w23 −w22w33)+
g4µ4(−w42w33w24 +w32w43w24 +w42
w23w34−w22w43w34−w32w23w44+w22w33w44))))−
g1µ1(v2(v3(v4w11 +g4µ4(w41w14 −w11w44))+
g3µ3(v4(w31w13 −w11w33)+
g4µ4(−w41w33w14 +w31w43w14 +w41w13w34−
w11w43w34 −w31w13w44 +w11w33 w44)))+
g2µ2(v3(v4(w21w12 −w11w22)+
g4µ4(−w41w22w14+w21w42 w14+w41w12w24−
w11w42w24 −w21w12w44 +w11w22w44))+
g3µ3(v4(−w31w22w13+w21 w32w13+w31w12w23−
w11 w32w23 −w21w12w33 +w11 w22w33)+
g4µ4(−w21w42w33w14 +w21w32w43w14+
w11w42w33w24 −w11w32w43w24+
w21w42w13w34 −w11w42w23w34−
w21w12w43w34 +w11w22w43w34+
w41(−w32w23w14 +w22w33w14 +w32w13w24−
w12w33w24 −w22w13w34 +w12w23w34)−
w21w32w13w44 +w11w32w23w44+
w21 w12w33w44 −w11w22w33w44+
w31(w42w23w14 −w22w43w14 −w42w13w24+
w12w43w24 +w22w13w44 −w12w23w44))))).
Such formulas are considered in paper [16].
2.2 Logistic function
The logistic function or logistic curve f(z) =
1
1+e−µ(z−θ)[7]. The sigmoid logistic func-
tion was introduced in a series of three pa-
pers by Pierre Francois Verhulst between
1838 and 1847, who devised it as a model
of population growth by adjusting the expo-
nential growth model. The sigmoid function
has the the characteristic properties:
1. monotonically increasing from zero to
unity;
2. possessing a unique inflection point,
[12].
-4-2 2 4
0.2
0.4
0.6
0.8
1.0
Figure 2: The sigmoid logistic function.
2.3 Critical points
The four-dimensional system has 4 eigenval-
ues.
•4D node. All eigenvalues are real and
have the same sign. The node is sta-
ble (unstable) when the eigenvalues are
negative (positive).
•4D star. All eigenvalues are equal.
The 4D star is stable (unstable) when
the eigenvalues are negative (positive).
•Saddle. All eigenvalues are real and
at least one of them is positive and at
least one is negative. Saddles are al-
ways unstable.
•Focus −Node. It has two real eigen-
values and a pair of complex-conjugate
eigenvalues, and all eigenvalues have
real parts of the same sign. The critical
point is stable (unstable) when the sign
is negative (positive).
•Node −Focus. It has two real nega-
tive eigenvalues and a pair of complex-
conjugate eigenvalues with positive real
part. The critical point is unstable.
•Saddle −Focus. Two real eigenval-
ues have different signs and complex-
conjugate eigenvalues with positive or
negative real part. The critical point is
unstable.
•Focus −Focus. Two pairs of
complex-conjugate eigenvalues. The
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.15
Inna Samuilik
E-ISSN: 2415-1521
115
Volume 10, 2022
critical point is stable when the signs
of real parts are negative. The critical
point is unstable when there is at least
one positive real part.
3 Materials and methods
Our consideration is geometrical. The main
intent is to use the 2D and 3D projections of
the 4D trajectories on different subspaces,
to construct the graphs of solutions for
understanding and managing the system.
Visualizations where possible, are provided.
The dynamics of Lyapunov exponents
are shown. Computations are performed
using Wolfram Mathematics, [12]. In the
article for Lyapunov exponents calculation
the package “lce.m for Mathematica” was
used, [13]. Another Wolfram Mathematica
program Lynch-DSAM.nb was also used
to check the correctness of Lyapunov
exponents calculation, [14], [15].
3.1 The example of four-
dimensional system
Consider the system (2) with the regulatory
matrix
W=
1220
−2 1 0 0
−2 0 0.7 2
0 0 −2 1
(5)
and µ1=µ3=µ4= 5, µ2= 15, v1=
v2=v3=v4= 1. In paper [9] θi, where
i= 1,2,3,4 are calculated as
θ1=w11 +w12 +w13 +w14
2,
θ2=w21 +w22 +w23 +w24
2,
θ3=w31 +w32 +w33 +w34
2,
θ4=w41 +w42 +w43 +w44
2.
θ1=2.5, θ2=−0.5, θ3= 0.35, θ4=−0.5.
The initial conditions are
x1(0) = 0.2; x2(0) = 0.15;
x3(0) = 0.3; x4(0) = 0.35.(6)
The critical point is (0.5; 0.5; 0.5; 0.5).
The characteristic equation for the
critical point is (4), where A=
−3.125, B = 32.2813, M =−39.8359
and L= 125.174. Solving the equation
we have λ1,2= 0.606091 ±2.15656iand
λ3,4= 0.956409 ±4.90201i. The type of
the critical point is unstable focus-focus.
The projection of 4D trajectories on two-
dimensional subspace (x1, x2) is in the
figure below.
0.0
0.2
0.4
0.6
0.8
x1
0.0
0.2
0.4
0.6
0.8
x2
Figure 3: The projection of 4D trajectories to
2D subspace (x1, x2).
The graphs of periodic solutions
(x1(t), x2(t), x3(t), x4(t)), the graphs of
periodic solutions (x1(t), x4(t)) and the
graphs of periodic solutions (x1(t), x4(t)) of
the system (2) with the regulatory matrix
(5) are shown in Figure 4, Figure 5 and
Figure 6.
The projection of 4D trajectories to 3D
subspace (x1, x2, x4) is shown in Figure 7.
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.15
Inna Samuilik
E-ISSN: 2415-1521
116
Volume 10, 2022
60
70
80
90
100
t
0.40
0.45
0.50
0.55
0.60
0.65
x1 x2x3x4
Figure 4: The graphs of periodic solutions
(x1(t), x2(t), x3(t), x4(t)) of the system (2) with
the regulatory matrix (5).
20
40
60
80
100
t
0.3
0.4
0.5
0.6
0.7
0.8
x1x2
Figure 5: The graphs of periodic solutions
(x1(t), x2(t)) of the system (2) with the reg-
ulatory matrix (5).
The dynamics of Lyapunov expo-
nents are shown in Figure 8. Lyapunov
exponents are LE1= 0.001, LE2=
−0.124, LE3=−0.127, LE4=−0.639 it
means (0,−,−,−). The system (2) with
the regulatory matrix (5) has periodic
solutions.
4 Conclusions
In this paper, we concern with mathema-
tical models of genetic networks. We have
considered the four-dimensional system.
Such a system has four characteristic
numbers which can be found from the
characteristic equation. Formulas for
20
40
60
80
100
t
0.3
0.4
0.5
0.6
0.7
x1x4
Figure 6: The graphs of periodic solutions
(x1(t), x4(t)) of the system (2) with the reg-
ulatory matrix (5).
0.0
0.5
1.0
X1
0.0
0.5
1.0
X2
0.0
0.5
1.0
X4
Figure 7: The projection of 4D trajectories to
3D subspace (x1, x2, x4).
the characteristic equation coefficients
(A, B, M, L) are considered. The projec-
tions of 4D trajectories to 2D subspace
and 3D subspace are considered using the
software Mathematica Wolfram. Also, the
graphs of periodic solutions of the system
(2) with the regulatory matrix (5) are
considered using the same software. As
we see the attractor can exist in the form
of an attracting closed trajectory (limit
cycle). In the future, we can perturbation
the regulatory matrix coefficients to have
other types of solutions for the system (2),
for example, chaotic solutions.
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.15
Inna Samuilik
E-ISSN: 2415-1521
117
Volume 10, 2022
0 500 1000 1500 2000 2500 3000
-1.0
-0.5
0.0
0.5
Steps
LCEs
Figure 8: LE1= 0.001, LE2=−0.124, LE3=
−0.127, LE4=−0.639.
References
[1] A. Lyapunov, G. Bagrinovskaja. On
Methodological Issues in Mathematical
Biology, Mathematical modeling in bi-
ology. - M, pp.5-18,1974.
[2] G.Demidenko, N. Kolchanov, V. Li-
hoshvai, J. Matushkin, S. Fadeev.
Mathematical modeling of regular con-
tours of gene networks, Journal of
computational mathematics and math-
ematical physics, vol.12, 22762295,
2004.
[3] B. D. MacArthur, P. S. Stumpf, R.
O.C.Oreffo. From mathematical mod-
eling and machine learning to clinical
reality, Principles of Tissue Engineer-
ing (Fifth Edition), pp.37-51, 2020.
[4] C. Furusawa, K. Kaneko. A
generic mechanism for adaptive
growth rate regulation. PLoS
Comput Biol 4(1), 2008: e3.
doi:10.1371/journal.pcbi.0040003
[5] H.D. Jong. Modeling and Simulation
of Genetic Regulatory Systems:
A Literature Review, J. Com-
put Biol. 2002;9(1):67-103, DOI:
10.1089/10665270252833208
[6] F. Sadyrbaev, I. Samuilik, V. Sen-
gileyev. On Modelling of Genetic Regu-
latory Networks. WSEAS Transactions
on Electronics, vol. 12, pp. 73-80, 2021
[7] I. Samuilik, F. Sadyrbaev, D.
Ogorelova. Mathematical modeling of
three-dimensional genetic regulatory
networks using logistic and Gompertz
functions, WSEAS Transactions on
systems and control, pp.101-107, 2022.
DOI: 10.37394/23203.2022.17.12
[8] J. Berro. Essentially, all models are
wrong, but some are useful-a cross-
disciplinary agenda for building use-
ful models in cell biology and bio-
physics, Biophysical Review, 10(6), pp.
16371647, 2018. DOI:10.1007/s12551-
018-0478-4.
[9] O. Kozlovska, F. Sadyrbaev. Models
of genetic networks with given prop-
erties, WSEAS Transactions on com-
puter reserch, pp. 43-49, 2022. DOI:
10.37394/232018.2022.10.6
[10] Y.Koizumi et al. Adaptive Virtual Net-
work Topology Control Based on At-
tractor Selection. J. of Lightwave Tech-
nology, (ISSN :0733- 8724), Vol.28
(06/2010), Issue 11, pp. 1720-1731
DOI:10.1109/JLT.2010.2048412.
[11] A.Das, A.B.Roy, Pritha Das. Chaos
in a three dimensional neural net-
work. Applied Mathematical Mod-
elling, 24(2000), 511-522.
[12] I. Samuilik, F. Sadyrbaev, S. Atslega.
Mathematical modelling of nonlinear
dynamic systems, Engineering for Ru-
ral Development, 21, pp. 172-178,2022.
[13] M. Sandri. Numerical calculation of
Lyapunov exponents, The Mathemat-
ica Journal, 1996.
[14] S. Lynch. Dynamical Systems with
Applications Using Mathematica.
Springer, 2017.
[15] I.Samuilik. Genetic engineering-
construction of a network of four
dimensions with a chaotic attractor,
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.15
Inna Samuilik
E-ISSN: 2415-1521
118
Volume 10, 2022
Vibroengineering Procedia, 44, pp.
66-70, 2022.
[16] F. Sadyrbaev, I. Samuilik. Remark on
four dimensional system arising in ap-
plications. Proceedings of IMCS of Uni-
versity of Latvia, 20(1), 2020. ISSN
1691 8134
Contribution of individual
authors to the creation of
a scientific article (ghost-
writing policy)
Author Contributions:
All authors have contributed equally to cre-
ation of this article.
Sources of funding for re-
search presented in a sci-
entific article or scientific
article itself
ESF Project No.8.2.2.0/20/I/003 Strength-
ening of Professional Competence of Dau-
gavpils University Academic Personnel of
Strategic Specialization Branches 3rd Call
Creative Commons Attri-
bution License 4.0 (Attri-
bution 4.0 International ,
CC BY 4.0)
This article is published under the terms of
the Creative Commons Attribution License
4.0
creativecommons.org/licenses/by/4.0/deed.en US
creativecommons.org/licenses/by/4.0/deed.en US
WSEAS TRANSACTIONS on COMPUTER RESEARCH
DOI: 10.37394/232018.2022.10.15
Inna Samuilik
E-ISSN: 2415-1521
119
Volume 10, 2022
