On trajectories of a system modeling evolution of genetic networks

1.   Introduction

Gene regulatory networks (GRN) exist in any cell of any living organism. Genetic networks are responsible for morphogenesis and adaptation to changes in the environment. Elements of these networks interact with each other, developing a general reaction of an organism to changes in the external environment. Biologists have created huge arrays of experimental data, which need to be analyzed, classified and utilized. For these purposes, mathematical models have been created. They are different, using mathematical theories such as Boolean algebras, graph theory, stochastic analysis and more. To follow the evolution of GRN in time, dynamical systems are best suited. Models of GRN, based on differential equations, can be treated qualitatively and computationally. On the basis of this treatment, predictions of the behaviors of GRN can be made. The structure of GRN and principles of functioning can be clarified by studying the phase spaces of dynamical systems in a model. Management and control over GRN can be made possible in some practically important cases. Genetic networks are important in the treatment of various diseases, such as multiple sclerosis, chemotherapy for gastric cancer and leukemia. GRN can be modeled mathematically, and in case of success, the model is to be studied, using analytical and computational approaches. Evolution in time can be traced, following the trajectories of a dynamical system, which describe the behavior of GRN. For this, knowledge of attracting sets is important. This system consists of several ordinary differential equations (ODE), containing multiple parameters. The interrelation between genes is described by the $n\times n$ regulatory matrix W. This system is not a particular one. Its modification ^[1], used in the theory of artificial neural networks, is known to be able to approximate to arbitrary accuracy any dynamical system ^[2].

We consider the n-dimensional dynamical system

where ${\mu }_{i} > 0$ and ${\theta }_{i}$ are parameters, and ${w}_{ij}$ are elements of the $n\times n$ regulatory matrix

Such systems arise in the theory of genetic regulatory networks ^{[3,4,5,6,7,8]}, neural networks ^[9,10] and telecommunication networks ^[11].

These systems contain nonlinear terms, represented by the function $\frac{1}{1+{e}^{-\mu z}}$, and the linear part also. The nonlinear function $f\left(z\right)$ is sigmoidal, that is, it is smooth, monotonically increasing from zero to unity and has a single inflection point. The sigmoidal functions are many, and the Hill and the Gompertz functions can be used as well. The argument of the function $f\left(z\right)$ can be replaced by the linear combinations of multiple variables ${x}_{1}, \dots, {x}_{n}.$ Thesevariables are interpreted as genes, expressing proteins and thuscommunicating with each other, developing a common response of thegene network to external and internal influences. In the theory of telecommunication networks ${x}_{i}$ are treated as source-destination pairs, and the optimal configuration of lightpaths is the problem to be solved.

Our attention to these problems stems from the fact that they are a focus of research on the border of mathematics and biology, and the classical mathematical approach should be corrected and adapted to the special needs of practicioners.

Some authors ^[12,13] interpret this system as a model for a not so large real genomic network, which is functioning in an organismsuffering from leukemia (a kind of blood cancer). Following Wang et al. ^[12], the genetic network can be in various states, and some of these states are associated with the disease, while some are considered as "normal". From the mathematical point of view, the respective system of order 60 (only!) has attractive equilibria (critical points), and some of them are conditionally "bad", and some are "normal". The current state of an organism is associated with the state vector $X\left(t\right) = \left({x}_{1}\left(t\right), \dots, {x}_{n}\left(t\right)\right)$, which consists of solutions of the autonomous system (1) represented by trajectories, which are curves connecting the current state $X\left(t\right)$, which depends on the initial conditions $X\left(0\right)$, and the "final" state, which is treated as an "attractor".

2.   Two-dimensional systems

It is reasonable to start with a system of the form (1), which contains only two equations and corresponds to a two-element network. The system has the form

and it defines the vector field $\left({x}_{1}, {x}_{2}\right)$, which can bevisualized. Any solution of the system can be viewed as a trajectory $\left({x}_{1}\left(t\right), {x}_{2}\left(t\right)\right)$, which is usually a curve in the phase plane $\left({x}_{1}, {x}_{2}\right).$ We will refer to this system as a 2D-system. Consequently, a system of three equations will be called a 3D-system and so on. Even in the 2D-case, the system contains eight parameters, and therefore the solutions of the system can exhibit various behaviors. Signs and values of the elements ${w}_{ij}$ of the regulatory matrix reflect the character of interrelations between elements ${x}_{i}$. The structure of the phase portraits heavily depends on the nullclines, which are defined by the equations

The solutions of the system (4) are critical points, which are called in some contexts "equilibria". The local analysis of solutions of 2D-systems can be conducted by considering the linearized (around a critical point) system and computing some values, called the characteristic values ^[14]. The three types of behaviors are shown in Figure 1(a)–(c). Usually the trajectories of a2D-system are attracted by stable equilibria. Sometimes the attracting set is a limit cycle. More complicated structures can be attractors as well.

In the theory of boundary value problems for ordinary differential equations, the quasilinear equations were paid special attention by researchers. Fixed point theorems of functional analysis were used to prove the existence of solutions to several quasi-linear problems (e.g., ^[15,16]). The system (1) is quasi-linear, and its 2D version (3) can be written in a convenient form,

The traditional boundary conditions that can be associated with (5) are

Since ${f}_{1}$ and ${f}_{2}$ in (5) are continuous and bounded, the above problems are quasi-linear. It follows from the general theory of boundary value problems (BVP) for ordinary differential equations (ODE) ^[15,16] that the problems (8) and (9) are solvable for any $A$ and $B$, since the corresponding homogeneous problems

(and for ${x}_{2}\left(0\right) = 0, {x}_{1}\left(T\right) = 0$ also) have the trivial solutions only.

This is not the case for the problems (6) and (7), since the homogeneous problems

(and for ${x}_{2}\left(0\right) = 0, {x}_{2}\left(T\right) = 0$ also) have nontrivial solutions in the form ${x}_{1}\left(t\right)\equiv 0, {x}_{2}\left(t\right) = {e}^{-t}$.

As to the periodic problem (10), it always has a solution in the form of a constant $\left({{x}^{*}}_{1}, {{x}^{*}}_{2}\right)$. This follows from the geometricalfact that both nullclines in (4) are located in the strips $0 < {x}_{1} < 1$ and $0 < {x}_{2} < 1$, respectively, and therefore must intersect. The more interesting case is to have a non-constant solution, but additional restrictions are required for this.

The most influential factor on the behavior of solutions is the matrix $W$, which is called the regulatory matrix (thus the abbreviation GRN).

Consider three matrices

and the respective 2D-systems of the form (3). In the theory of genetic networks these types of interrelations between genes are called, respectively, activation, inhibition and mixed activation-inhibition. To be definite, let other parameters be ${\mu }_{1} = {\mu }_{2} = 5, {\theta }_{1} = \frac{{w}_{11}+{w}_{12}}{2}, {\theta }_{2} = \frac{{w}_{21}+{w}_{22}}{2}$. These choices of $\theta$ put a critical pointat the central location.

The cross-points of nullclines are marked by the small black disks. Only one critical point is in Figure 1(c). The phase portrait (on the larger scale) of the same system is depicted in Figure 1(d). It is clear that the problem (9) is not solvable for some values of $A, B$.

All the boundary conditions above mean that solutions are sought which connect subspaces of the extended phase space.

For the purposes of modeling genetic networks, other conditions fit better. Look at Figure 1(a). Both side critical points have their basins of attraction, that is, the sets of initial conditions which are connectable with the respective critical point. Both basins of attraction are separated by separate axes of the middle critical point, which is a saddle. At least in the unit quadrant ${Q}_{2},$ there are no trajectories connecting points from different basins of attraction. This means that the BVP with the boundary sets

is not solvable in ${Q}_{2}$. This boundary problem is over-determined, since the number of conditions (four) is greater than the number of equations (two).

2.1.   The alternative formulation of a boundary value problem

Let $A$ be an attracting set for the system (3) ^[17]. This means that there exists a vicinity $U$ of $A$ such that for any trajectory Γ with the initial conditions in $U$, this trajectory stays in $U$ eventually, and the distance d(Γ, $A$) tends to zero as t goes to infinity. The corresponding boundary condition can be written as $\left({x}_{1}\left(+{\mathit{\infty}} \right), {x}_{2}(+{\mathit{\infty}})\right)\in A.$ Of course, astable equilibrium and a limit cycle are attractors. Let us formulate the problem:

where $T$ is some number, which may be different for different trajectories. Since the second condition is looked at as depending on $T$, it can be written in a neutral form $\left({x}_{1}\left(t\right), {x}_{2}\left(t\right)\right)\in A.$

Example 1.

Let the 2D-system be with the matrix ${W}_{a}$ as in (13). Let ${A}_{1}$ be the upper stable equilibrium (it is a stable node). The problem is

This problem is solvable for any $P\left(0\right)$ belonging to the basin ofattraction of the point ${A}_{1}$. The problem has no solutions if $P\left(0\right)$ is in the basin of attraction of the lower stable equilibrium ${A}_{2}$.

Example 2.

Consider the 2D-system corresponding to the matrix ${W}_{m}$ in (13). Let $L$ be the limit cycle which is depicted in Figure 1(c). The problem is

3.   Three-dimensional systems

The classical boundary value problems for the third-order systems of the form (1) can be

where

The functions ${f}_{i}$ are sigmoidal functions of the form as in (1). They are continuous and bounded, and in accordance with the general theory ^[15,16], they are solvable, if the linear homogeneous problem has only the trivial solution.

For instance, the problems (18), (21) have a solution for any $A$, $B$, $C$ since the homogeneous problem has only the trivial solution for any T.

The periodic problem is a special one. It has a solution since any system of the form (1) has at least one critical point (the constant solution).

As in the case of 2D systems, the alternative boundary value problem can be posed in the form

Example 3.

Consider the system (1), where ${\mu }_{1} = {\mu }_{3} = 5, {\mu }_{2} = 15, {\theta }_{1} = {\theta }_{3} = 1.5, {\theta }_{2} = -0.25$, and the regulatory matrix is

This system has exactly three critical points, which are calculated and depicted as the cross-points of the nullclines in Figure 2. The first one is $\left(0.001;0.305;0.001\right).$ The linearization around it in a standard way yields the linear system. Let the 3 by 3 coefficient matrix of the linearized system be denoted as M. The characteristic equation det(M−λE) = 0 can be written in the form

where $A = -6.18;B = -9.33;C = -4.16$. Solving the equation, we have ${\lambda }_{1} = -4.18;{\lambda }_{2} = -1.003;{\lambda }_{3} = -0.9916$. The type of this critical point is a stable node.

The characteristic equation for the critical point $\left(0.5;0.5;0.5\right)$ is (29), where $A = -4.25;B = 8.56;C = 29.39$. Solving the equation, we have ${\lambda }_{1} = -4.75;{\lambda }_{2} = -2.25;{\lambda }_{3} = 2.75$. The type of this critical point is a saddle.

The characteristic equation for the critical point $\left(0.99;0.69;0.99\right)$ is (28), where $A = -6.18;B = -9.33;C = -4.16$. Solving the equation, we have ${\lambda }_{1} = -4.18;{\lambda }_{2} = -1.003;{\lambda }_{3} = -0.9916$. The type of this critical point is a stable node.

Example 4.

Consider ${\mu }_{1} = {\mu }_{3} = 5, {\mu }_{2} = 15, {\theta }_{1} = 1.5;{\theta }_{2} = 1, {\theta }_{3} = -0.5$. The regulatory matrix of the system (1) is

There are exactly three critical points. The characteristic equation for the critical point $\left(0.5;0.5;0.5\right)$ is (28), where $A = 3.25;B = -7.69;C = 17.36$. Solving the equation, we have ${\lambda }_{1} = 2.75;{\lambda }_{\mathrm{2, 3}} = 0.25\pm 2.5i$. The type of the critical point is a 3D unstable focus-node. The nullclines and the graph of solutions are depicted in Figure 3. The self-excited attractor and the dynamics of Lyapunov exponents are depicted in Figure 4.

Example 5.

Consider the system

where ${v}_{1} = 0.65, {v}_{2} = 0.42, {v}_{3} = 0.1, {\mu }_{1} = {\mu }_{2} = 7, {\mu }_{3} = 13, {\theta }_{1} = 0.5, {\theta }_{2} = 0.3, {\theta }_{3} = 0.7$. It is a slight modification of the system from the works ^[19,20]. The regulatory matrix of the system (31) is

The initial conditions are

There is one critical point. The characteristic equation for the critical point $\left(0.37;1.59;0.22\right)$ is (29), where $A = -1.16;B = -0.42;C = -0.69$. Solving the equation, we have ${\lambda }_{1} = -1.26;{\lambda }_{\mathrm{2, 3}} = 0.05\pm 0.74i$. The critical point isof the 3D unstable saddle-focus type. The system is chaotic in the sense that solutions exhibit non-regular behavior.

4.   Conclusions

The boundary value problems of a new kind arose from studies of mathematical models associated with genetic, neural and telecommunication networks. The behavior of GRN on bounded time intervals can be described by systems of ordinary differential equations of the form (1). The interrelation between elements of GRN is encrypted in the regulatory matrix W. The individual properties of genes are described by the parameters µ and θ. In the mathematical model (1), the unit cube Q is an invariant set, and trajectories cannot exit it. The nullclines can (and they do) intersect only in Q. The equilibria (critical points) exist. It is possible that they are attractive, and cases where all the equilibria are non-attractive are possible also. Then, more complicated attractive sets exist, and more often they are periodic ones. The boundary value problems of the forms (16), (17), (27) arise naturally in this context. The boundary value problems on finite intervals also can be studied in light of suggestions in the paper ^[12]. The problem of redirecting a trajectory from the basin of attraction of a "bad" attracting set to a "normal" one, changing the adjustable parameters, is of great practical importance.

Acknowledgments

This work is supported by the following project: ESF Project No. 8.2.2.0/20/I/ "Strengthening of Professional Competence of Daugavpils University Academic Personnel of Strategic Specialization Branches 3rd Call."

Conflicts of interest

The authors declare no conflict of interest.