Modeling Networks of Four Elements

Abstract

In this article, fourth-order systems of ordinary differential equations are studied. These systems are of a special form, which is used in modeling gene regulatory networks. The nonlinear part depends on the regulatory matrix W, which describes the interrelation between network elements. The behavior of solutions heavily depends on this matrix and other parameters. We research the evolution of trajectories. Two approaches are employed for this. The first approach combines a fourth-order system of two two-dimensional systems and then introduces specific perturbations. This results in a system with periodic attractors that may exhibit sensitive dependence on initial conditions. The second approach involves extending a previously identified system with chaotic solution behavior to a fourth-order system. By skillfully scanning multiple parameters, this method can produce four-dimensional chaotic systems.

1. Introduction

Gene networks in living organisms regulate basic processes and participate in the control of vital functions. They are extensively studied by teams of investigators of various kinds: biologists, chemists, physicists, geneticists, statisticians, and mathematicians. Mathematical models are actively used to study gene networks. Models of gene regulatory networks are diverse. They can be deterministic or stochastic. In deterministic models, the future states of a system are determined by the initial state and external input. Due to the possible chaotic behavior of solutions, this does not make it easier to analyze future states. On the other hand, in real-world systems, the stochastic effects may play an important role. Stochastic models are not always necessary; it depends on the modeled system. Generally, genetic networks possess mechanisms that contribute to their overall stability. Maintaining a balance between complexity and stability is important for properly functioning biological systems. When the switch-like behavior among groups of genes is observed, discrete models (such as Boolean algebra and graph theory-based) may be preferable. The evolution of a network in time can be modeled using ordinary differential equations. Details on this topic are available in review articles [1,2,3,4,5].

Our focus is on models formulated as dynamical systems. There are numerous models. Some are mentined in review article [6]. The difference model is described where the increase in the expression level of gene i at time t + Δt is dependent on the sum of weighted expression levels of other genes at time t. An artificial neural network (ANN) method is mentioned that is capable of describing the dynamic behavior of genetic regulatory networks (GRNs). An inhibitory model of three ordinary differential equations (sometimes referred to as the “Goodwin Oscillator”) was analyzed in [7]. Typically, ODE models of gene networks incorporate nonlinearities that are expressed using sigmoidal functions. In [8], an algorithm for qualitative simulation of gene regulatory networks with steep sigmoidal functions is provided. In [9], there is a focus on “a mechanistic understanding of the gene regulation underneath”. For this, a single-cell model is studied. The author claims that “resolving gene expression in individual cells can allow subtle behaviors to be detected, and by achieving these feats with thousands to millions of data points per experiment, single-cell approaches promise to properly grasp the complexity of biological systems”. In [10], GRNs are treated as chemical reaction networks. For this, the ODE apparatus is needed. For example, the kinetics of the lac operon interpreted as a CRN is described by a system of eleven ODEs. In [11], “Several commonly used network modeling paradigms are surveyed with emphasis on their practical use in systems biology research”. When speaking about ODEs modeling genetic networks, the authors say that “the central idea is to model the rate of change in each reaction as a set of differential equations. These differential equations are coupled, which means the rate of change in the expression of a gene is affected by the current expression level of other genes. When formulating an ODE model, it is important to choose the right balance between model complexity, identifiability, and interpretability”. It is also noted that “it is usually impractical to model all the physical biomolecular processes in a GRN (such as transcription, translation, translocation, and so on). Many have proposed to model functional relationships among genes instead”.

The concept of neural ODEs that belongs to the field of deep learning and its application for modeling GRNs are discussed in [12]. The modeling framework based on neural ODEs is proposed and tested in a series of experiments involving realistic GRNs. In [13] the authors present “an extended model for gene regulatory networks with time delayed negative feedback, which is described by delay differential equations”. In [14] “the finite volume discretization scheme for partial integro-differential equations (PIDEs) describing the temporal evolution of protein distribution in gene regulatory networks is proposed”. This ODE-based model “allows the application of the theory of nonnegative and compartmental systems for the qualitative analysis of the approximating dynamics”. The goal of article [15] “is to develop genetic networks with temporal delays, which are crucial for genetic regulation because slow biochemical processes like gene transcription and translation need time to occur”.

System (2) in the article under review has a long history. In 1972, H. Wilson and J. Cowan studied the interaction of two populations of neurons, active and inhibitory [16]. The qualitative analysis of the phase plane for a two-dimensional system was performed. A simplified system in the form (1) was studied qualitatively and numerically in [17,18]. The structure of the phase plane and the behavior of solutions were in focus. The 3D system was considered in several sources providing new results on the possible behavior of solutions. An interesting question of the role of stochasticity in cellular reaction dynamics resulting in the selection of a stable gene expression pattern was investigated in [19]. Gene regulation with positive and negative feedback in an n-dimensional model (2) was studied in [18]. System (2) was treated as “a neural network, where each node represents a particular gene and the wiring between the nodes defines regulatory interactions”.

The second group of related results, focusing on the construction of systems with a block diagonal regulatory matrix W and considering systems in search of chaos, is associated with the work [20] and the literature therein. In [21], the 6D example of system (2) was examined.

The third group of problems concerns biomedicine-oriented applications of systems of the form (2). In [22,23], system (2) was considered as a model suitable for considering and finding ways of treating leukemia (a kind of blood cancer). In a model, the problem is interpreted as going the system state trajectories to a “wrong” attractor. The ways of redirecting trajectories to a “normal” attractor were discussed in [24,25].

The fourth group of problems involves identifying chaotic attractors within system (2). Although numerous chaotic attractors may exist in systems of type (2), their discoveries are rare. References [26,27] document the same 3D chaotic attractor. The chaotic attractors in a 3D system of type (2) are discussed in [28,29].

The dynamics and evolution of networks can be captured by n-dimensional ODE-based models. Attracting sets play a crucial role in these models, as they define the configuration of phase spaces and enable important conclusions about trajectories and the properties of the modeled networks. Each network element is represented by a time-dependent value of x_i(t). The components of the network are connected by certain relations, which are described by the regulatory matrix W, which is built into the system. The current state of the network is represented by the vector X(t) = (x₁(t),…, x_n(t)), which contains all the elements of the network. The structure of the proposed system ensures that the system’s trajectories, represented by the vectors X(t), are confined within a specific bounded region of the phase space. Consequently, the system can have attractors, which are sets that attract trajectories. Describing attractors in a system of differential equations is crucial for forecasting the network’s future states. Available examples suggest that attractors can take the form of critical points, periodic attractors, and more complex sets emerging from the chaotic behavior of the system’s solutions. Various types of attractors have been discovered in two-dimensional models, three-dimensional models, and higher-dimensional individual systems. The study of attractors is not an easy task, since GRN systems contain many parameters, the number of which increases significantly as the dimensionality of the system increases. In this paper, a four-dimensional (4D for short) model of a gene network containing 28 parameters is studied. The work consists of two parts, in which the 4D model is studied using the approaches previously used by the authors. In the first part, systems are built, which are initially a combination of two 2D systems with known properties. These two 2D systems are made dependent, combining into one 4D system, with the help of special techniques. Initially, two two-dimensional systems with periodic attractors are merged into a four-dimensional system. The same process is then applied to another pair of 2D systems, each possessing a periodic attractor. Next, both pairs are included in a family of 4D systems, depending on the parameter p, which varies in the interval [0, 1]. For p = 0 we have the first 4D system; for p = 1 the second one is obtained. Scanning the resulting family of 4D systems by changing the parameter p, we obtained an attractor, the projection of which on the first three axes is depicted in Figure 1. In the next example, the two 2D systems are similarly combined into one 4D system. The first 2D system has an attractor in the form of a limit cycle with clockwise rotation. The second 2D system has an attractor as a limit cycle with counterclockwise rotation. The resulting 4D system is perturbed by adding nonzero elements to the regulatory matrix. An attractor appears in the resulting 4D system. The sensitive dependence of solutions on the initial data is observed by constructing the Lyapunov exponents [30]. More examples of that kind are provided. The second part of the article deals with the same subject but uses another approach. Previously obtained and published results on the 3D systems are used [28,29]. By adding suitable components to the 3D matrices used in the works, a 4D system is obtained, which has a visualizable 4D attractor. The behavior of solutions is chaotic. The Lyapunov curve analysis is provided. Useful information on the interpretation of the Lyapunov spectrum is in [30]. Information on the best understanding of gene regulatory networks and chaotic behaviors of solutions can be found in [31]. The properties of nonlinear systems and the role of attractors are covered in detail in [32].

Despite being one of the most extensively studied areas in higher mathematics, the field of ordinary differential equations continues to expand significantly. This growth is driven by its widespread practical applications and its importance in education. The latter stems from the fact that linear and nonlinear ODE theory offers numerous instances of analytical problem solving, including the theory of linear ODE with constant coefficients.

Typically, ordinary differential equations (ODEs) cannot be solved analytically and require computational methods. To illustrate how attractors can emerge, we present two approaches, both utilizing numerical tools.

The second section covers preliminaries and basic knowledge of systems in the forms (1), (5), and (10). The third section examines four-dimensional systems with attractors that may exhibit sensitive dependence on initial conditions. Section 4 provides new examples of 4D systems exhibiting chaotic solutions behavior.

2. Preliminaries

First, we consider the two-dimensional (2D) systems of ordinary differential equations (ODE) of the form

(1)$\left\{\begin{array}{c}{x^{\prime }}_1={\displaystyle \frac{1}{1+e^{-{\mu }_1\left(w_{11}{x}_1+w_{12}{x}_2-{\theta }_1\right)}}}-v_1{x}_1,\\ {x^{\prime }}_2={\displaystyle \frac{1}{1+e^{-{\mu }_2\left(w_{21}{x}_1+w_{22}{x}_2-{\theta }_2\right).}}}-v_2{x}_{2.}\end{array}\right.$

This system is the simplification [17,18] of the system that first appeared in the research by Wilson–Cowan [16] on the dynamics of interactions between populations of excitatory and inhibitory neurons.

The right side contains two terms, a nonlinear one and a linear term. Without the nonlinear term, the system possesses only exponentially decreasing solutions, and it is not interesting (but physically explainable: no evolution without interactions).

Generally, μ_i are positive numbers, greater than 1, θ_i are adjustable parameters, and v_i are positive, but in most cases they are set to 1. The parameters of the GRN have the following biological interpretations:

v_i—degradation of the i-th gene expression product; w_ij—the connection weight or strength of control of gene j on gene i. Positive values of w_ij signify activating influences, whereas negative values denote repressing influences; θ_i—the impact of external stimuli on gene i is reflected in its ability to modulate the gene’s responsiveness to activating or repressing factors [26].

An invariant set in a dynamical system is a set of points in the system’s phase space that, once entered, the system will remain within for all future time steps. The invariant set for system (1) is the rectangle $Q_2=\left\{\left(x_1, x_2\right): 0<x_1<{\displaystyle \frac{1}{v_1}}, 0<x_2<{\displaystyle \frac{1}{v_{12}}} \right\}.$ Further, the n-dimensional system

(2)$\left\{\begin{array}{c}{x^{\prime }}_1={\displaystyle \frac{1}{1+e^{-{\mu }_1({\mathrm{w}}_{11}{x}_1+\dots +w_{1n}{x}_n-{\theta }_1)}}}-v_1{x}_1,\\ \dots \dots \dots \\ {x^{\prime }}_n={\displaystyle \frac{1}{1+e^{-{\mu }_n({\mathrm{w}}_{n1}{x}_1+\dots +w_{nn}{x}_n-{\theta }_n)}}}-v_{\mathrm{n}}{x}_{\mathrm{n}},\end{array}\right.$

was used extensively in the mathematical modeling of gene regulatory networks and telecommunication networks [33].

The basic facts about systems (1) and 3D version of system (2) can be found in [24] and references therein, and in [28], theorems 1 and 2.

Here are some of the facts:

(1). The above-defined set $Q_2$ is invariant;

(2). In this set there is at least one critical point;

(3). The nullclines, defined by the equations

(3)$\left\{\begin{array}{c}0={\displaystyle \frac{1}{1+e^{-{\mu }_1\left(w_{11}{x}_1+w_{12}{x}_2-{\theta }_1\right)}}}-v_1{x}_1,\\ 0={\displaystyle \frac{1}{1+e^{-{\mu }_2\left(w_{21}{x}_1+w_{22}{x}_2-{\theta }_2\right).}}}-v_2{x}_{2,}\end{array}\right.$

intersect only in the invariant set; therefore, all possible critical points are located in the invariant set.

Multiple examples of the study of the critical points of system (2) can be found in [24,28,29].

The behavior of trajectories of system (1) essentially depends on the coefficients w_ij. There are three main modes of behavior, namely, activation, inhibition, and mixed case (which results in a rotating vector field) in system (1). They are characterized, respectively, by the following matrices:

(4)$W_{active}=\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right), W_{inhib}=\left(\begin{array}{cc}1& -1\\ -1& 1\end{array}\right), W_{inhib}=\left(\begin{array}{cc}1& 1\\ -1& 1\end{array}\right).$

3. Four-Dimensional Systems of Special Structures

3.1. Regulatory Matrices

Consider the four-dimensional system of the form

(5)$\left\{\begin{array}{c}{x^{\prime }}_1=f_1({\mathrm{w}}_{11}{x}_1+\dots +{\mathrm{w}}_{14}{x}_4-{\theta }_1)-v_1{x}_1,\\ {x^{\prime }}_2=f_2({\mathrm{w}}_{21}{\mathrm{x}}_1+\dots +{\mathrm{w}}_{24}{\mathrm{x}}_4-{\theta }_2)-v_2{x}_2,\\ {x^{\prime }}_3=f_3({\mathrm{w}}_{31}{\mathrm{x}}_1+\dots +{\mathrm{w}}_{34}{\mathrm{x}}_4-{\theta }_3)-v_3{x}_3,\\ {x^{\prime }}_4=f_4({\mathrm{w}}_{41}{\mathrm{x}}_1+\dots +{\mathrm{w}}_{44}{\mathrm{x}}_4-{\theta }_4)-v_4{x}_4,\end{array}\right.$

where ${\mathit{f}}_{\mathit{i}}\left(\mathit{z}\right)$ are sigmoidal, that is, they are monotonically increasing from 0 to 1 and have only one point of inflection. The matrix of coefficients

(6)$W=\left(\begin{array}{cccc}{w}_{11}& w_{12}& w_{13}& w_{14}\\ w_{21}& w_{22}& w_{23}& w_{24}\\ w_{31}& w_{32}& w_{33}& w_{34}\\ w_{41}& w_{42}& w_{43}& w_{44}\end{array}\right)$

is called the regulatory matrix. Let W be of the form

(7)$W_{41}=\left(\begin{array}{cccc}{w}_{11}& w_{12}& 0& 0\\ w_{21}& w_{22}& 0& 0\\ 0& 0& w_{33}& w_{34}\\ 0& 0& w_{43}& w_{44}\end{array}\right),$

where the non-zero elements are chosen so that the two 2D systems each have periodic attractors. This is always possible. The resulting 4D attractor is the product of two 2D attractors. Its 3D projections onto subspaces can be observed.

Let

(8)$W_{42}=\left(\begin{array}{cccc}0& 0& w_{13}& w_{14}\\ 0& 0& w_{23}& w_{24}\\ w_{31}& w_{32}& 0& 0\\ w_{41}& w_{42}& 0& 0\end{array}\right)$

be another matrix with a block structure, where any 2D system corresponding to a block has a limit cycle as an attractor.

Consider the resulting matrix

(9)$W_{res}=p{W}_{41}+\left(1-p\right)W_{42},$

where p ϵ [0, 1]. This matrix is W₄₁ for p = 1, and it is W₄₂ for p = 0. System (5) is expected to have an attractor.

Choose the sigmoidal function f(z) = 1/(1 + exp(−μ(w z − θ))) and write system (5) as

(10)$\left\{\begin{array}{c}{x^{\prime }}_1={\displaystyle \frac{1}{1+e^{-{\mu }_1({\mathrm{w}}_{11}{x}_1+\dots +w_{14}{x}_4-{\theta }_1)}}}-v_1{x}_1,\\ {x^{\prime }}_2={\displaystyle \frac{1}{1+e^{-{\mu }_2({\mathrm{w}}_{21}{x}_1+\dots +w_{24}{x}_4-{\theta }_2)}}}-v_2{x}_2,\\ {x^{\prime }}_3={\displaystyle \frac{1}{1+e^{-{\mu }_3({\mathrm{w}}_{31}{x}_1+\dots +w_{34}{x}_4-{\theta }_3)}}}-v_3{x}_3,\\ {x^{\prime }}_4={\displaystyle \frac{1}{1+e^{-{\mu }_4({\mathrm{w}}_{41}{x}_1+\dots +w_{44}{x}_4-{\theta }_4)}}}-v_4{x}_4.\end{array}\right.$

Let matrix (7) be

(11)$W_{41}=\left(\begin{array}{cccc}0.6& 1& 0& 0\\ -1& 0.6& 0& 0\\ 0& 0& 0.6& 2\\ 0& 0& -2& 0.6\end{array}\right)$

It consists of two 2D blocks, and these blocks correspond to a 2D system of the form (1) where all four ${\mu }_i=8$. These systems exhibit attractive periodic solutions of periods T₁ ≈ 3.78 and T₂ ≈ 1.88, respectively, as measured on solution graphs obtained in Wolfram Mathematica.

Similarly, construct matrix (8) as

(12)$W_{42}=\left(\begin{array}{cccc}0& 0& 1& 1\\ 0& 0& -1& 1\\ 0.548& 1& 0& 0\\ -1& 0.548& 0& 0\end{array}\right)$

The non-zero 2D blocks correspond to two 2D systems that exhibit stable periodic solutions, with periods of approximately T₃ » 6.82 and T₄ » 3.41, respectively.

3.2. A Parameterized Family of 4D Systems

Construct the 4D matrix following the scheme in (9).

It is

(13)$W_{big}=\left(\begin{array}{cccc}0.6p& p& (1-p)& (1-p)\\ -p& 0.6p& -(1-p)& (1-p)\\ 0.548(1-p)& (1-p)& 0.6p& 2p\\ -(1-p)& 0.548(1-p)& -2p& 0.6p\end{array}\right)$

Let p ϵ [0, 1]. For p = 0.95, we obtain the 4D attractor in system (10). The 3D projections of the attractor can be visualized using Mathematica software. The 3D projection on the xyz-subspace (it is called for convenience the 123 projection) is depicted in Figure 1a. The Lyapunov curves are presented in Figure 1b.

Now, let us construct the 4D matrices and examine the examples. We will create several attractors, which are perturbations of the periodic attractors for the unperturbed matrices.

3.3. Counter-Rotating Systems

Set the parameters μ_i = 4 (for the remainder of this section). Let the regulatory matrix (counter-rotating) be defined as follows:

(14)$W_{CR}=\left(\begin{array}{cccc}2& 1& 0& 0\\ -1& 2& 0& 0\\ 0& 0& 2& -1\\ 0& 0& 1& 2\end{array}\right)$

The 2D block in the upper left corner corresponds to the system

(15)$\left\{\begin{array}{c}{x^{\prime }}_1={\displaystyle \frac{1}{1+e^{-4\left(2x_1+x_2-1.5\right)}}}-x_1,\\ {x^{\prime }}_2={\displaystyle \frac{1}{1+e^{-4\left(-x_1+2x_2-0.5\right).}}}-x_{2.}\end{array}\right.$

It has the limit cycle depicted in Figure 2a. The parameters θ₁ and θ₂ are set to 1.5 and 0.5, respectively, to keep a critical point (unstable focus) at the center.

Similarly, the 2D block in the right lower corner corresponds to the system

(16)$\left\{\begin{array}{c}{x^{\prime }}_3={\displaystyle \frac{1}{1+e^{-4\left(2x_3-x_4-0.5\right)}}}-x_3,\\ {x^{\prime }}_4={\displaystyle \frac{1}{1+e^{-4\left(x_3+2x_4-1.5\right).}}}-x_{4.}\end{array}\right.$

which has the limit cycle depicted in Figure 2b. The parameters θ₃ and θ₄ are set to 0.5 and 1.5, respectively, to keep a critical point (unstable focus) at the center.

The vector field for one system rotates in the opposite direction to the vector field for the other system.

Both limit cycles when combined into one 4D system

(17)$\left\{\begin{array}{c}{x^{\prime }}_1={\displaystyle \frac{1}{1+e^{-4\left(2x_1+x_2-1.5\right)}}}-x_1,\\ {x^{\prime }}_2={\displaystyle \frac{1}{1+e^{-4\left(-x_1+2x_2-0.5\right)}}}-x_2,\\ {x^{\prime }}_3={\displaystyle \frac{1}{1+e^{-4\left(2x_3-x_4-0.5\right)}}}-x_3,\\ {x^{\prime }}_4={\displaystyle \frac{1}{1+e^{-4\left(x_3+2x_4-1.5\right)}}}-x_{4.}.\end{array}\right.$

with matrix (14), produce the 4D attractor. The 3D projections of this attractor are visualized in Figure 3a. These projections were obtained by visualization of some trajectories of system (17). The trajectory with initial conditions (0, 0, 0, 0) converges to the black component of the graph in Figure 3a. The trajectory with initial conditions (0, 1, 0, 0) reaches the blue oval in Figure 3a, while the trajectory starting from the point (0, 0, 1, 0) ends up at the red oval, which is the limiting set.

Consider the perturbed system

(18)$\left\{\begin{array}{c}{x^{\prime }}_1={\displaystyle \frac{1}{1+e^{-4\left(2x_1+x_2+0.5 x_4-1.75\right)}}}-x_1,\\ {x^{\prime }}_2={\displaystyle \frac{1}{1+e^{-4\left(-x_1+2x_2-0.5\right)}}}-x_2,\\ {x^{\prime }}_3={\displaystyle \frac{1}{1+e^{-4\left(2x_3-x_4-0.5\right)}}}-x_3,\\ {x^{\prime }}_4={\displaystyle \frac{1}{1+e^{-4\left({0.5 x_1+x}_3+2x_4-1.75\right)}}}-x_{4.}.\end{array}\right.$

with the regulatory (counter-rotating perturbed) matrix

(19)$W_{CRP}=\left(\begin{array}{cccc}2& 1& 0& 0.5\\ -1& 2& 0& 0\\ 0& 0& 2& -1\\ 0.5& 0& 1& 2\end{array}\right)$

Several trajectories of system (18), tending to a 4D attractor, are depicted in Figure 3b. The results of testing system (18) for sensitive dependence on initial conditions, as measured by the Lyapunov exponents, are displayed in Figure 3c.

The positive Lyapunov curve in Figure 3c refers to the sensitive dependence of solutions of system (16) on the initial data. In the two figures below the solutions and trajectories emanating from almost the same initial point are depicted. Further, the trajectories diverge significantly.

The solutions are sensitive, dependent on the initial data. It can be seen in Figure 4a that at the initial stage, oscillation of the red solution is more extensive than that of the blue one.

Remark 1.

In our examples, the parameters θ_i are chosen as the half-sum of the elements in the respective row of the related regulatory matrix. This will keep the critical point at the central location (0.5, 0.5, 0.5, 0.5). For instance,

(20)$\left\{\begin{array}{c}0={\displaystyle \frac{1}{1+e^{-4\left(2 0.5+0.5+0.5 0.5-1.75\right)}}}-0.5,\\ 0={\displaystyle \frac{1}{1+e^{-4\left(-0.5+2 0.5-0.5\right)}}}-0.5,\\ 0={\displaystyle \frac{1}{1+e^{-4\left(2 0.5-0.5-0.5\right)}}}-0.5,\\ 0={\displaystyle \frac{1}{1+e^{-4\left(0.5 0.5+0.5+2 0.5-1.75\right)}}}-0.5.\end{array}\right.$

for system (18). More on this can be found in [34].

3.4. Different Rates of Activation

Let the regulatory matrix be

(21)$W_{CRk}=\left(\begin{array}{cccc}k& 1& 0& 0\\ -1& k& 0& 0\\ 0& 0& 1.2& -1\\ 0& 0& 1& 1.2\end{array}\right)$

Setting θ₁ = 0.5(k + 1), θ₂ = 0.5(−1 + k), θ₃ = 0.5(2 − 1), θ₄ = 0.5(1 + 2), we keep the critical point at (0.5, 0.5, 0.5, 0.5). The standard linearization at (0.5, 0.5) shows that the 2D system corresponding to the block

(22)$\left(\begin{array}{cc}k& 1\\ -1& k\end{array}\right)$

has a critical point with zero real part and non-zero imaginary one. The value k = 1 in (21) is a bifurcation point for the Andronov—Hopf bifurcation. Compose the 4D system with the regulatory matrix

(23)$W_{CR1}=\left(\begin{array}{cccc}1& 1& 0& 0\\ -1& 1& 0& 0\\ 0& 0& 1.2& -1\\ 0& 0& 1& 1.2\end{array}\right)$

The activation coefficients on the main diagonal are different, respectively, 1 and 1.2. The 3D projection of the respective 4D attractor along with several trajectories is shown in Figure 5.

(24)$W_{CR1P}=\left(\begin{array}{cccc}1& 1& 0& 0\\ -1& 1& 0& 0\\ 0& 0& 1.2& -1\\ 0.5& 0& 1& 1.2\end{array}\right)$

Adding an element of 0.5 to the matrix (23) leads to the sensitive dependence in system (10) with the previously specified parameters and regulatory matrix (24).

4. Four-Dimensional Systems Based on 3D Chaotic Systems

In this section, we present 4D systems exhibiting chaotic behavior. Previously, we constructed several 3D chaotic systems. Various experiments were conducted to identify suitable values for the matrix elements that result in chaotic solutions. We form 4D systems from chaotic 3D ones by adding additional matrix elements: w₁₄, w₂₄, w₃₄, w₄₄, and parameters: v₄, μ₄, θ₄. Example 1 is based on the 3D chaotic attractor found and described in [28]. All key characteristics, including the analysis of critical points, the location and shape of nullclines, and Lyapunov exponents, were obtained.

4.1. Example 1

The 4D system was formed using 3D system that was studied in [28]. This 3D system was shown to have three unstable critical points and a chaotic attractor. A comparison with the Chua circuit was made.

Consider the resulting 4-dimensional system (10) and the regulatory matrix

(25)$W_{ex1}=\left(\begin{array}{cccc}0.06& -0.825& 0.19& 0\\ 1.915& 1.76& -0.32& 0\\ 3.7& 0.128& 0.985& 0\\ 1& 1& 0& 1\end{array}\right)$

v₁ = 1, v₂ = 1, v₃ = 1, v₄ = 1; μ₁ = 4; μ₂ = 5.02, μ₃ = 8, μ₄ = 8; θ₁ = (w₁₁ + w₁₂ + w₁₃ + w₁₄)/2, θ₂ = (w₂₁ + w₂₂ + w₂₃ + w₂₄)/2, θ₃ = (w₃₁ + w₃₂ + w₃₃ + w₃₄)/2, θ₄ = (w₄₁ + w₄₂ + w₄₃ + w₄₄)/2. Let initial conditions be (0.592; 0.85; 2; 1). The solutions for system (10) with matrix (25) are depicted in Figure 6. It can be observed that the solutions exhibit irregular behavior. For a nonlinear dynamic system, the number of Lyapunov exponents coincides with the dimension of the phase space. There are four such indicators in our example. The results of testing system (10) on the Lyapunov exponents for the initial data (0.592; 0.85; 2; 1) are (0.05; 0.00; −0.45; −0.78), where LE_1 > 0, LE_2 = 0, LE_3 < 0, and LE_4 < 0. A positive Lyapunov exponent indicates chaotic behavior in the system [30].

The projections of this attractor onto the three-dimensional subspaces are shown in Figure 7.

4.2. Example 2

We consider the 4-dimensional system obtained from the 3-dimensional one in [29] by adding elements for the fourth column and the fourth row of the regulatory matrix W. The 3D system studied in [29] was constructed using the methodology by L. Shilnikov. This methodology requires the existence of saddle-focus points satisfying some additional conditions. As a result, rich dynamical behavior of solutions can be observed, including chaotic behavior. The history, references, and construction details can be found in [29].

Consider the 4-dimensional system (10) with the regulatory matrix

(26)$W_{ex2}=\left(\begin{array}{cccc}0& 0.262& -7.12& 0.1\\ 1.46& 0& 4& 0\\ 0.14251& -1& 2& 0\\ 1& 0& 0& 1\end{array}\right)$

and v₁ = 0.708, v₂ = 0.307, v₃ = 0.767, v₄ = 1; μ₁ = 5.63, μ₂ = 4.538, μ₃ = 4, μ₄ = 1; θ₁ = −4, θ₂ = 4.44, θ₃ = 0.5, θ₄ = 0.5. The initial conditions are (0.4; 0.9; 0.4; 0).

The Lyapunov exponents are (0.04; 0.00, −0.64; −1.11) with the initial condition (0.4; 0.9; 0.4; 0), where LE₁ > 0, LE₂ = 0, LE₃ < 0, and LE₄ < 0. The solutions of system (10) with regulatory matrix (26) are shown in Figure 8. The three-dimensional projections of this attractor are shown in Figure 9.

5. Discussion

The study of gene networks is one of the important tasks of biomathematics. Systems of differential equations make it possible to model networks in motion and monitor their evolution. In this case, the problem becomes purely mathematical, but its tasks are corrected with the help of the biological data that are obtained. The system studied in this article is an interesting mathematical object. It is quasi-linear, i.e., it is the sum of linear and nonlinear parts. This property, together with the range of values of the nonlinearity (usually a sigmoidal function with the range (0, 1)), ensures the presence of an invariant set in the phase space and, as a result, the existence of an attractor(s). Attractors determine the future states of the system (and network). Attractors can be regular, the properties of which are well traced, and chaotic.

For this study, it is useful to construct examples of attractors for systems in question, where possible. In the first part of this work, such attractors are constructed as geometric objects obtained by combining systems of lower dimensions. In this case, well-readable systems of dimension 4 are obtained. This dimension was chosen for the study because it is the beginning of the non-visualizable field of study of gene systems. Several attractors have been obtained, the structure of which can be imagined by analyzing three-dimensional projections on the coordinate axes.

Attractors whose properties are close to chaotic were obtained, demonstrating the sensitivity of solutions to the initial data. An attractor has been received, the projections of which look like the union of closed curves in a 3-dimensional space. The tendency of the trajectories of the main 4-dimensional system to such an attractor depends on the initial data. Solutions tend to a selected attractor (based on the initial data). The biological counterpart of this behavior, tending to a selected periodic regime, should be studied and interpreted. These attractors are well visualized. In the second part of the work, two examples of chaotic attractors for 4-dimensional systems are constructed. For this, previously studied 3-dimensional systems with 3-dimensional chaotic attractors are used. The challenge is to augment these 3D systems to 4D systems in a nontrivial way, increasing the regulatory matrix W by one row and one column and saving the property of having a chaotic attractor. It is to be mentioned that other (concerning the elements of matrix W) parameters of the system are also involved in the process, for example, coefficients vi at the linear terms of a system.

There are some additional comments on the ongoing text. Examples of 3D systems that have a chaotic attractor were constructed in [28,29]. Paper [28] uses the idea of the Chua circuit, which suggests that the system has 3 critical points (equilibria) of a certain type. The Chua circuit is a relatively simple-looking three-dimensional (3D) system that is almost linear with only one piecewise linear function on the right. The construction of a similarly organized 3D system of GRN-type (Genetic Regulatory System type) required work-consuming efforts. In [29], the 3D chaotic system of the GRN type was constructed following the scheme proposed by Shilnikov. This scheme required the existence of a saddle-focus type equilibrium (critical point) with a certain balance between the characteristic numbers. The extension of these two 3D examples to the 4D examples with chaotic attractors was performed by adding a few elements to the matrices W, making it four-dimensional, and checking the result.

We suggest that our research advances the field in several directions. First, we have considered the parameter-dependent family of 4D systems, choosing afterward the appropriate system with the needed behavior. Second, we have considered 4D systems combined with two 2D systems with definite behavior of solutions. These two 2D systems were asked to have periodic solutions (limit cycles) with certain properties. In the second part of the work, the chaotic attractors were constructed for the 4D systems. This is a remarkable fact since even 3D examples of chaotic attractors for GRN-type systems are rare (for algebraic systems, the examples are many).

In our examples, two interacting subsystems, each representing a periodic process, are combined into a larger system. For high-dimensional systems where subsystems can be identified, a similar approach may approximate the larger system, potentially leading to more efficient analysis.

There are biological guidelines for selecting parameters in system (2). Altering μ impacts a gene’s individual property to accept inputs. Experimental data may offer insights into the type of interaction within a network or subsystems. Therefore, the matrix W, or an approximation, can be reconstructed. In a model, particularly for the 2D case, the activatory matrix W (where all elements are positive) generally results in stable equilibria. The same is true for the inhibition, but equilibria are located differently. In the mixed case, where there is a balance between activation and inhibition, stable periodic attractors may exist.

Periodic attractors were obtained in examples of the first part. They presumably may correspond to the repeating and rhythmic processes in real networks.

Lyapunov exponents relate to the chaotic behavior of solutions, or, at least, to sensitive dependence on the initial data. This goes beyond bistability and/or oscillatory behaviors. The bistability, or switching, of the behavior of genes has a natural mathematical counterpart, namely, certain bifurcations of solutions under the change in parameters. Mathematically, this can occur when the trajectory of a dynamical system crosses the border between basins of attraction of two neighboring attractors. The oscillations are a focus of our approach. Several kinds of oscillatory behaviors can be met even in 4D systems. Oscillations in GRN-type systems are very common. They are stimulated by a balance between activatory and inhibitory populations of genes.

Generally, GRN-type systems are stable (rough, structurally stable). Therefore, the influence of noise is not extensive and can be even neglected in many studies of GRN. The influence of unknown parameters can be hidden in the parameter θ. It can significantly affect the main properties of a system. On the other hand, this allows for management and control of the network, which is important for applications.

6. Conclusions

Periodic attractors constructed in 2D systems can be used to build 4D attractors in 4D systems. By perturbation of the 4D regulatory matrices, multiple 4D systems can be obtained exhibiting interesting properties. The multi-fragmented attractors can be constructed for 4D systems with the property that solutions of the system with a perturbed regulatory matrix tend to a fragment of the 4D attractor. The choice of the fragment depends on the initial values of a solution. On the other hand, 4D systems with the desired behavior of solutions can be constructed by adding one dimension (one suitably selected equation) to a three-dimensional system with previously detected remarkable properties. In this way, 4D systems with chaotic behavior can be obtained. Several examples of this kind are provided, utilizing the results from papers [28,29]. Moreover, the two chaotic systems obtained in this paper can be used to construct systems of higher dimensions, such as eight dimensions and beyond. However, the main challenge remains interpreting the obtained behaviors in the context of real biological networks.

Author Contributions

Conceptualization, F.S. and O.K.; methodology, F.S. and O.K.; formal analysis, writing—F.S. Section 1, Section 2 and Section 3, O.K. Section 4. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No new data was obtained during this research.

Acknowledgments

The authors thank both reviewers for pointing out possible directions for future research.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1 (a) The 123 projection of the attractor corresponding to matrix (13), p = 0.95; (b) Lyapunov curves.

Figure 2 The vector fields for (a) system (15) and (b) system (16).

Figure 3 (a) The 123 projections of the trajectories of system (17). Black is for the trajectory with the start point at (0, 0, 0, 0), blue is for the initial conditions (0, 1, 0, 0), red is for the initial conditions (1, 0, 0, 0), and dashed for (0, 0, 1, 0); (b) the 123 projections of the trajectories of system (18) with the initial conditions (0, 0, 0, 0) (black), (0, 1, 0, 0) (blue), (1, 0, 0, 0) (red), (0, 0, 1, 0) (dashed); (c) the Lyapunov curves for the trajectory of system (18). The initial conditions are (0.4, 0.6, 0.4, 0.6). The Lyapunov exponents are {{0.253554, −0.0246819, −0.556984, −0.671894}, 2.41091}. The last number is the Kaplan–York–Lyapunov fractal dimension.

Figure 4 (a) The ${\mathit{x}}_1\left(\mathit{t}\right)$ components of two solutions of system (18) with the initial conditions (0.4, 0.6, 0.4, 0.6) (red) and (0.41, 0.6, 0.4, 0.6) (blue); (b) The 123 projections of two trajectories of system (18) with the initial conditions (0.4, 0.6, 0.4, 0.6) (red) and (0.41, 0.6, 0.4, 0.6) (blue).

Figure 5 (a) The 123 projections of several trajectories, going to the attractor, matrix (24); (b) Lyapunov curves for the trajectory starting at (0.4, 0.6, 0.4, 0.6), matrix (23), Lyapunov exponents (0.0118, −0.0167, −0.0385, −0.2905); (c) Lyapunov curves for the trajectory starting at (0.4, 0.6, 0.4, 0.6), matrix (24), Lyapunov exponents (0.2535, −0.0247, −0.5569, −0.6719).

Figure 6 The solutions of system (8) with regulatory matrix (25).

Figure 7 The 3D projections of the 4D attractor in system (10) with matrix (25): (a) the 234 projection of the 4D attractor; (b) the 124 projection of the 4D attractor; (c) the 134 projection of the 4D attractor; and (d) the 123 projection of the 4D attractor.

Figure 8 The solutions of system (10) with regulatory matrix (26).

Figure 9 The 3D projections of the 4D attractor in system (10) with matrix (26): (a) the 234 projection of the 4D attractor; (b) the 124 projection of the 4D attractor; (c) the 123 projection of the 4D attractor; and (d) the 134 projection of the 4D attractor.

References

1. Alakwaa, F.M. Modeling of Gene Regulatory Networks: A Literature Review. J. Comput. Syst. Biol.; 2015; 1, 1. [DOI: https://dx.doi.org/10.15744/2455-7625.1.102]

2. Yang, B.; Chen, Y. Overview of Gene Regulatory Network Inference Based on Differential Equation Models. Curr. Protein Pept. Sci.; 2020; 21, pp. 1054-1059. [DOI: https://dx.doi.org/10.2174/1389203721666200213103350] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32053072]

3. Bhattacharya, P.; Raman, K.; Tangirala, A.K. Discovering adaptation-capable biological network structures using control-theoretic approaches. PLOS Comput. Biol.; 2022; 18, e1009769. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1009769] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35061660]

4. Emmert-Streib, F.; Dehmer, M.; Haibe-Kains, B. Gene regulatory networks and their applications: Understanding biological and medical problems in terms of networks. Front. Cell Dev. Biol.; 2014; 2, 38. [DOI: https://dx.doi.org/10.3389/fcell.2014.00038]

5. McCallum, G.; Potvin-Trottier, L. Using Models to (Re-)Design Synthetic Circuits. Synthetic Gene Circuits; Methods in Molecular Biology; Menolascina, F. Humana: New York, NY, USA, 2021; Volume 2229, [DOI: https://dx.doi.org/10.1007/978-1-0716-1032-9_3]

6. Vijesh, N.; Chakrabarti, S.K.; Sreekumar, J. Modeling of gene regulatory networks: A review. J. Biomed. Sci. Eng.; 2013; 6, pp. 223-231. [DOI: https://dx.doi.org/10.4236/jbise.2013.62A027]

7. Smith, H. A Classical Gene Regulatory Network. Available online: https://sites.science.oregonstate.edu/~deleenhp/teaching/spring06/MAP4484/hal-cyclicgene.pdf (accessed on 12 March 2025).

8. Ironi, L.; Panzeri, L.; Plahte, E. An Algorithm for Qualitative Simulation of Gene Regulatory Networks with Steep Sigmoidal Response Functions. Algebraic Biology; Horimoto, K.; Regensburger, G.; Rosenkranz, M.; Yoshida, H. Lecture Notes in Computer Science Springer: Berlin/Heidelberg, Germany, 2008; Volume 5147, [DOI: https://dx.doi.org/10.1007/978-3-540-85101-1_9]

9. Maizels, R.J. A dynamical perspective: Moving towards mechanism in single-cell transcriptomics. Philos. Trans. R. Soc. B Biol. Sci.; 2024; 379, 20230049. [DOI: https://dx.doi.org/10.1098/rstb.2023.0049]

10. Barbuti, R.; Gori, R.; Milazzo, P.; Nasti, L. A survey of gene regulatory networks modelling methods: From differential equations, to Boolean and qualitative bioinspired models. J. Membr. Comput.; 2020; 2, pp. 207-226. [DOI: https://dx.doi.org/10.1007/s41965-020-00046-y]

11. Ho, J.W.K.; Charleston, M.A. Network modelling of gene regulation. Biophys. Rev.; 2010; 3, pp. 1-13. [DOI: https://dx.doi.org/10.1007/s12551-010-0041-4]

12. Hossain, I.; Fanfani, V.; Fischer, J.; Quackenbush, J.; Burkholz, R. Biologically informed NeuralODEs for genome-wide regulatory dynamics. Genome Biol.; 2024; 25, pp. 1-29. [DOI: https://dx.doi.org/10.1186/s13059-024-03264-0]

13. Öztürk, D.; Atay, F.M.; Özbay, H. Chaos in gene regulatory networks: Effects of time delays and interaction structure. Chaos Interdiscip. J. Nonlinear Sci.; 2024; 34, 033102. [DOI: https://dx.doi.org/10.1063/5.0172767]

14. Vághy, M.A.; Otero-Muras, I.; Pájaro, M.; Szederkényi, G. A Kinetic Finite Volume Discretization of the Multidimensional PIDE Model for Gene Regulatory Networks. Bull. Math. Biol.; 2024; 86, pp. 1-26. [DOI: https://dx.doi.org/10.1007/s11538-023-01251-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38253903]

15. Cao, Y.; Chandrasekar, A.; Radhika, T.; Vijayakumar, V. Input-to-state stability of stochastic Markovian jump genetic regulatory networks. Math. Comput. Simul.; 2023; 222, pp. 174-187. [DOI: https://dx.doi.org/10.1016/j.matcom.2023.08.007]

16. Wilson, H.R.; Cowan, J.D. Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons. Biol. Cybern.; 1972; 13, pp. 1-18. [DOI: https://dx.doi.org/10.1016/S0006-3495(72)86068-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/4332108]

17. Noonburg, V.W. Differential Equations: From Calculus to Dynamical Systems; 2nd ed. AMS/MAA TEXTBOOKS: Providence, RI, USA, 2019; Volume 43.

18. Vohradsky, J. Neural network model of gene expression. FASEB J.; 2001; 15, pp. 846-854. [DOI: https://dx.doi.org/10.1096/fj.00-0361com]

19. Furusawa, C.; Kaneko, K. A generic mechanism for adaptive growth rate regulation. PLOS Comput. Biol.; 2008; 4, e3. [DOI: https://dx.doi.org/10.1371/journal.pcbi.0040003]

20. Samuilik, I. Genetic engineering—Construction of a network of four dimensions with a chaotic attractor. Vibroengineering PROCEDIA; 2022; 44, pp. 66-70. [DOI: https://dx.doi.org/10.21595/vp.2022.22829]

21. Kopp, M.; Samuilik, I. A New 6D Two-wing Hyperchaotic System: Dynamical Analysis, Circuit Design, and Sinchronization. Chaos Theory Appl.; 2024; 6, pp. 273-283. [DOI: https://dx.doi.org/10.51537/chaos.1513080]

22. Wang, L.-Z.; Su, R.-Q.; Huang, Z.-G.; Wang, X.; Wang, W.-X.; Grebogi, C.; Lai, Y.-C. A geometrical approach to control and controllability of nonlinear dynamical networks. Nat. Commun.; 2016; 7, 11323. [DOI: https://dx.doi.org/10.1038/ncomms11323]

23. Cornelius, S.P.; Kath, W.L.; Motter, A.E. Realistic control of network dynamics. Nat. Commun.; 2013; 4, pp. 1-9. [DOI: https://dx.doi.org/10.1038/ncomms2939]

24. Ogorelova, D.; Sadyrbaev, F.; Sengileyev, V. Control in Inhibitory Genetic Regulatory Network Models. Contemp. Math.; 2020; 1, pp. 393-400. [DOI: https://dx.doi.org/10.37256/cm.152020538]

25. Ogorelova, D.; Sadyrbaev, F.; Samuilik, I. On Targeted Control over Trajectories of Dynamical Systems Arising in Models of Complex Networks. Mathematics; 2023; 11, 2206. [DOI: https://dx.doi.org/10.3390/math11092206]

26. Das, P.; Das, A.; Roy, A.B. Chaos in a three-dimensional general model of neural network. Int. J. Bifurc. Chaos Appl. Sci. Eng.; 2002; 12, pp. 2271-2281. [DOI: https://dx.doi.org/10.1142/S0218127402005820]

27. Mukherjee, S.; Palit, S.K.; Bhattacharya, D.K. Is one dimensional Poincaré map sufficient to describe the chaotic dynamics of a three dimensional system?. Appl. Math. Comput.; 2013; 219, pp. 11056-11064. [DOI: https://dx.doi.org/10.1016/j.amc.2013.04.043]

28. Kozlovska, O.; Sadyrbaev, F.; Samuilik, I. A New 3D Chaotic Attractor in Gene Regulatory Network. Mathematics; 2023; 12, 100. [DOI: https://dx.doi.org/10.3390/math12010100]

29. Kozlovska, O.; Sadyrbaev, F. In Search of Chaos in Genetic Systems. Chaos Theory Appl.; 2024; 6, pp. 13-18. [DOI: https://dx.doi.org/10.51537/chaos.1380419]

30. Sandri, M. Numerical Calculation of Lyapunov exponents. Math. J.; 1996; 6, pp. 78-84.

31. Sayama, H. Introduction to the Modeling and Analysis of Complex Systems; Open SUNY: Albany, NY, USA, 2015.

32. Perko, L. Differential Equations and Dynamical Systems; 3rd ed. Springer: New York, NY, USA, 2001.

33. Koizumi, Y.; Miyamura, T.; Arakawa, S.; Oki, E.; Shiomoto, K.; Murata, M. Adaptive virtual network topology control based on attractor selection. J. Light. Technol.; 2010; 28, pp. 1720-1731. [DOI: https://dx.doi.org/10.1109/JLT.2010.2048412]

34. Kozlovska, O.; Sadyrbaev, F. Models of genetic networks with given properties. WSEAS Trans. Comput. Res.; 2022; 10, pp. 43-49. [DOI: https://dx.doi.org/10.37394/232018.2022.10.6]