Remarks on GRN-type systems

Systems of ordinary differential equations that appear in gene regulatory networks theory are considered. We are focused on asymptotical behavior of solutions. There are stable critical points as well as attractive periodic solutions in two-dimensional and three-dimensional systems. Instead of considering multiple parameters (10 in a two-dimensional system) we focus on typical behaviors of nullclines. Conclusions about possible attractors are made.


Introduction
We consider systems of ordinary differential equations (ODE) of the form, x 0 1 ¼ 1 1 þ e Àl 1 ðw 11 x 1 þw 12 x 2 þ:::þw 1n xnÀh 1 Þ À c 1 x 1 ; x 0 2 ¼ 1 1 þ e Àl 2 ðw 21 x 1 þw 22 x 2 þ:::þw 2n xnÀh 2 Þ À c 2 x 2 ; ::: x 0 n ¼ 1 1þe Àl n ðw n1 x 1 þw n2 x 2 þ:::þwnn xnÀhn Þ À c n x n ; where all parameters w ij and h i are real numbers, l i and c i are positive numbers. The function f ðzÞ ¼ 1 1þe Àlz is a sigmoidal function, 1 which is defined for all z 2 R, strictly increasing from zero to unity and has exactly one inflection point. Some other sigmoidal function, such as the Gompertz function 2 f ðzÞ ¼ e Àe Àlz , also can be used in the above system [1]. The system (1) is often used in models of various real-world phenomena. It has appeared (for n = 2) in the paper by Wilson and Cowan [2] as a simplified model of interaction of two populations of neurons. The n dimensional version is believed to describe the interrelations in genetic regulatory networks (GRN in short), that is, probably, the most popular object related to the system (1). A GRN can be thought of as a complicated network existing in living cells that instantly remodel themselves. Every element of this network, gene, interacts with others and the character of this interaction is described by the so-called regulatory matrix W. The positive entries w ij of W mean activation of x i by x j . The negative ones mean inhibition. The zero entries are for no relation. The dynamics of the entire GRN is described by the system (1). Future states of a network are of primary interest. To understand them, one has to investigate possible attractors of the system (1). There are various mathematical models of GRN, for instance, models, formulated in terms of graph theory, Boolean algebras, etc. The possibility of studying future network states is a significant advantage of models based on systems of differential equations. For a description of genes and genome, we invite the reader to consult papers and book chapters in [3][4][5][6].
One of the main features of GRN is the ability to adapt to changes in the environment. Understanding of the mechanism of adaptation is important for designers of artificial networks. The telecommunication networks are of this kind. In the works by several authors, for instance, [7,8], the self-organization of GRN is considered as a source of ideas for designers of practical networks.
It was proposed [8] to use the attractor selection as the principal mechanism for the creation of virtual network topology (VNT), which could realize the rapid response of telecommunication networks to environmental changes without human participation. When creating the attractor selection mechanism for network resource management, at first the so called regulatory matrix W should be defined which shows relationships between node pairs, that is, how each node pair affects each other including itself. Likely as in GRN networks, three types of influence exist, namely, activation, inhibition and no relation, corresponding to positive, negative and zero entries w ij of the regulatory matrix.
Both mentioned networks can be huge. Therefore, the respective mathematical models, formulated in terms of differential equations, also include systems of high dimensionality. The standard mathematical (qualitative and quantitative) analysis of such systems, except particular cases, seems to be difficult if it is standard. The creation of new concepts and ideas for solving this challenging problem is an actual task for mathematicians. The identification, description, systematization of attractors for systems of differential equations of high dimensionality is, therefore, an urgent task.
The classical theory of Poincaré and Bendixson exist for the two-dimensional systems. Even this theory has a number of unanswered questions, especially concerning limit cycles as nontrivial attracting sets. The three-dimensional systems have richer possibilities for the existence of attractors, other than the stable equilibria. The description and full systematization of three-dimensional attractors seems to be far from completion. Systems of the form (1), besides of their practical importance in light of the above discussion, are a very convenient object for finding and investigation of attractors of unusual nature.
In what follows we made some remarks and discuss some observations, concerning the two-dimensional and threedimensional systems of the form (1). We provide the number and possible location of stable equilibria for typical particular cases. We consider the conditions for an equilibrium to be of the type focus. We discuss the process of bifurcation of parameters leading to the Andronov -Hopf bifurcation and emerging the isolated periodic solutions in two-dimensional and three-dimensional cases.

Boundary value problems
System (1) is a particular quasi-linear system of ordinary differential equations. These systems were intensively studied with respect to the relative boundary value problems before. The general theorem [9,10], adapted for the case of system (1) states, that the boundary value problem (1), where A i are n Â n matrices with real entries, t i are time moments in a definite interval, u is a continuous function of m variables, is solvable (classical continuously differentiable solution), if the function in the right hand side of (2) is bounded and the homogeneous problem, . . .
where x(t) = (x 1 (t), . . ., x n (t)), has only the trivial (zero) solution. This theorem is easy to apply for particular boundary conditions. Therefore the problem of existence of a solution, subject to some boundary conditions, not necessarily twopoint ones, can be treated effectively. The problems on infinite intervals are specific. They often require some knowledge about the attractors of a given system. Consider the two-dimensional system, together with the boundary conditions, where a, A, b, B are real numbers. The homogeneous problems, (6) and (9), (7) may have nontrivial solutions. Generally the problem (9), (6) (and/or (9), (7)) has not a solution. This is easy to see. Consider the system, together with the conditions, A solution x 1 (t) is uniquely defined by the initial condition x 1 (0) = A. Therefore a solution (x 1 (t), x 2 (t)) to the problem (10), (11) does not exist if B is not appropriate. The above reasoning is valid also for studying the problem with periodic or antiperiodic boundary conditions. For instance, the problem (13), is solvable for any positive T. This does not mean, however, that there always exists a closed orbit for the autonomous system, since constant solutions also are periodic, and due to the properties of a sigmoidal nonlinearities in (5) at list one critical point exists.

Two-dimensional cases
Consider the system, that contains 10 parameters. We wish to discuss the structure of attractors for this system. We avoid the standard analysis of critical points due to multiple parameters and instead consider few variants which depend on the geometry and location of nullclines. These features depend of course on parameters but explicitly the values of parameters are not involved. It seems, that this system was not studied in such context by standard tools of a phase plane and critical points. The few related references are [1,11]. The regularity matrix is, The equilibria (critical points) are to be found from the system, The type of a (non-degenerate) critical point ðx Ã 1 ; x Ã 2 Þ is to be detected in a standard way concerning the linearized system, u 0 1 ¼ Àc 1 u 1 þ l 1 w 11 g 1 u 1 þ l 1 w 12 g 1 u 2 ; u 0 2 ¼ Àc 2 u 2 þ l 2 w 21 g 2 u 1 þ l 2 w 22 g 2 u 2 ; ( ð16Þ where, Let us list the basic facts about the system (13). The nullclines of the system (13) are defined by the relations (15). Due to properties of sigmoidal functions the nullclines are located in the strips 0 < x 1 < 1 c 1 and 0 < x 2 < 1 c 2 respectively. Therefore all the critical points (that are cross-points of nullclines) lie in the rectangle domain Q ¼ 0; 1 There is at least one critical point. The inspection of the boundary of Q shows that the vector field ðx 0 1 ; x 0 2 Þ is directed inward. Any trajectory of the system (13) stays in Q eventually.
The characteristic equation for a critical point ðx Ã where, q ¼ c 1 c 2 À g 1 c 2 w 11 l 1 À g 2 c 1 w 22 l 2 À g 1 g 2 w 12 w 21 l 1 l 2 þ g 1 g 2 w 11 w 22 l 1 l 2 : Up to nine critical points are possible in the system (13) [1].

Graphical analysis
Even the two-dimensional system (13) contains 10 parameters and standard analysis of all possible cases is rather long. Therefore it is convenient to treat the system graphically using the nullclines method.
In the first case there is one equilibrium and the vector field ðx 0 1 ; x 0 2 Þ is directed towards the critical point which is a sink. The second case is visualized in Figure 2. The nullclines are crossed at the first equilibrium and have a tangent point at the second one. The upper (simple) critical point is a sink. It can be shown [1] that the lower critical point is degenerate in the sense that one of the characteristic numbers k is zero. It is attractive, however.
Due to S-shape of nullclines there are at most three equilibria. This is shown in Figure 3. The side equilibria are sinks and the middle one is a saddle point. If w 11 and w 22 are not zeros more critical points are possible [1].

Mixed case
Consider system (13) provided that c 1 = c 2 = 1, w 11 = w 22 = 0, and w 12 and w 21 have opposite signs. The nullclines together with the vector field are depicted in Figure 4. If w 12 is positive then the rotation of the vector field is in the counter clock-wise direction.
Proposition 3.1 There is a single critical point of the type stable focus for the system (13) if w 11 = w 22 = 0 and w 12 and w 21 have opposite signs. Proof in [1]. We get from (19) to (21) that the characteristic equation is, k 2 þ 2k þ ð1 À g 1 g 2 w 12 w 21 l 1 l 2 Þ ¼ 0: The roots are À1 AE ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi g 1 g 2 w 12 w 21 l 1 l 2 p : Since w 12 w 21 is negative, a unique critical point is focus. If w 12 is negative and w 21 positive then (Fig. 5) the rotation of vector field is in the clock-wise direction. E. Brokan and F. Sadyrbaev: 4open 2020, 3, 8 Remark 3.1. There is no a closed orbit in the system (13) if w 11 = w 22 = 0. The divergence of the vector field, defined by (13), is not zero and the Bendixson criterium ( [12], Ch. X, [5]) for closed orbits in the region Q is not satisfied. If w 11 and w 22 are allowed to be nonzero, the limit cycle is possible [1,11] (see Fig. 6). Imagine that w 11 = w 22 = k > 0 The limit cycle appears as a result of Andronov -Hopf bifurcation. The unique critical point then changes the stability to instability as k increases from zero to unity and the real parts of the complex conjugate characteristic numbers k(k) become positive passing through zero.

Three-dimensional cases
Consider the three-dimensional system, 11 x 1 þw 12 x 2 þw 13 xnÀh 1 Þ À c 1 x 1 ; x 0 2 ¼ 1 1 þ e Àl 2 ðw 21 x 1 þw 22 x 2 þw 23 x n Àh 2 Þ À c 2 x 2 ; The nullclines for the system are defined by the relations, The critical points for the system (23) are the cross points of the nullclines. The linearized system for any critical point where, The nullcline analysis can work well in the three-dimensional case also. Since a critical point is the cross-point of all nullclines, the number of equilibria can be established by combination of visual observations and checking by roots finding programs. For the case under consideration the following situation is of particular interest. There exists at least one critical point in systems (23). If there are several non-attractive equilibria (saddles or unstable focuses) then an unusual behavior of solutions in the bounded domain is expected. The famous Lorenz and Rössler systems are of this kind. The Lorenz system has three critical points (a saddle and two unstable focuses) for particular values of parameters. The Rössler system has one unstable focus and a saddle. If there is exactly one critical point of this kind, then the trajectories cannot be trapped by it.
We provide an example of one critical point that is unstable focus (that is, two of three characteristic numbers are complex conjugate with positive real part).

Conclusions
The nullcline method together with other tools can be effectively applied in the study of asymptotic behavior of solutions and the structure of attractors for systems of low dimensions. The conclusions made on the basis of nullcline analysis can generate the hypotheses and suggestions for higher dimensional cases. Some questions are unanswered up to now. We have examples of a system (23) with a single critical point of the type unstable focus (two complex conjugate characteristic numbers have positive real part). It is possible also to have exactly three critical points in such configuration: two unstable focuses and a saddle point.