Limit cycles in models of circular gene networks regulated by negative feedback loops

Background The regulatory feedback loops that present in structural and functional organization of molecular-genetic systems and the phenomenon of the regulatory signal delay, a time period between the moment of signal reception and its implementation, provide natural conditions for complicated dynamic regimes in these systems. The delay phenomenon at the intracellular level is a consequence of the matrix principle of data transmission, implemented through the rather complex processes of transcription and translation.However, the rules of the influence of system structure on system dynamics are not clearly understood. Knowledge of these rules is particularly important for construction of synthetic gene networks with predetermined properties. Results We study dynamical properties of models of simplest circular gene networks regulated by negative feedback mechanisms. We have shown existence and stability of oscillating trajectories (cycles) in these models. Two algorithms of construction and localization of these cycles have been proposed. For one of these models, we have solved an inverse problem of parameters identification. Conclusions The modeling results demonstrate that non-stationary dynamics in the models of circular gene networks with negative feedback loops is achieved by a high degree of non-linearity of the mechanism of the autorepressor influence on its own expression, by the presence of regulatory signal delay, the value of which must exceed a certain critical value, and transcription/translation should be initiated from a sufficiently strong promoter/Shine-Dalgarno site. We believe that the identified patterns are key elements of the oscillating construction design.

(Continued from previous page) Fig. 1 Simplest autorepressilators: a negatively autoregulated one-gene network, b its graphical scheme, c graph of two-genes network, d variety of sub-networks in a trigenic network, e graph of n-elements gene network stimulation), Lrp (leucine-responsive regulatory protein) etc [6][7][8][9]. A lot of similar examples is known for the eukaryotes as well. This type of regulation is characteristic for RB (retinoblastoma) gene [10], hes family basic helix-loop-helix factor (bHLH) transcription factors, which are important for coordinated somite segmentation [11,12], DNA binding homeobox transcription factor Nanog, which is involved in embryonic stem cell proliferation, renewal, and pluripotency [13,14], and many other genes [15][16][17][18], etc.
The next stage of regulatory circuits' complexity is represented by the two-gene systems, where each gene encodes expression regulator of another gene (Fig. 1c). Here, the examples of negative regulation circuit are pairs of regulatory factors Sox3-Snail and DDX21-Snail. They are repressors of each other expression and control formation of ecto-and mesoderm at early stage of embryo development, [19]. For the Bacillus subtilis, such properties have products of the genes ComK (competence transcription factor) and kre (ComK repressor) which participate in differentiation of cellular phenotypes [20]. There are similar systems regulated by RNA molecules, for example non-coding RNA MicF and global regulator of E.coli metabolism Lrp act negatively on expression of each other [21,22].
The following class of complexity is composed by trigenic regulatory circuits, where each gene encodes transcriptional factor, or another protein, and regulation is the circular one. An example of such a system is negative regulatory circuit AMPK-mTOR-NRF2 which acts in regulation processes of autophagy [23].
Studies of dynamical characteristics of circular gene networks, in particular networks regulated by negative feedback loops, are motivated not only by their prevalence in molecular genetic systems, which is important per se, but also by their role in dynamics of different biological processes. Even the simplest regulatory circuits incorporated into complicated gene networks can determine dynamics of the process in the large, as it was shown for the gene HES7, whose oscillating expression dynamics determines process of vertebrata somitogenesis and is related with its autorepression properties [24]. Now, we study only simplest circular regulatory circuits which can be recovered in natural gene networks. This is specified by the rapid increase of the amount of subsystems when the amount of genes grows. For example, if the number n of gene elements in a circular gene network equals 3, then the number s(3) of its subnetworks equals 13 (Fig. 1d). If n = 4, then s(4) equals 201, if n = 5, then s(5) = 9390; similarly s(6) = 1531336, s(7) = 880492496, s(8) = 1792477159408 etc. So, the mathematical analysis of these subsystems seems to be technically unfeasible.
We do not consider here cyclic networks regulated by mixed feedback mechanisms, which are positive and negative at the same time. The models of these gene networks possess not only a high oscillatory potential [25][26][27] but also the possibility of complex, chaotic dynamics formation [28][29][30][31], including the hyperchaotic one [32][33][34]. Here we consider mathematical models of natural circular gene networks. Our studies are based both, on numerical simulations of these networks and on analytic mathematical methods. We describe limit cycles in simplest circular gene networks models and elaborate algorithms of construction of these cycles.

Analysis of limit cycles in a model of the simplest autorepressilator
The model of simplest circular regulatory circuit which acts by the negative feedback loop is represented by an equation Description and interpretation of this equation and that of positive parameters α and h are given below in the "Methods" section. It was shown in [35] that the Eq. (1) has exactly one positive equilibrium point If τ = 0, then for all values of parameters α and h this equilibrium is stable, and all trajectories of the Eq. (1) tend to x 0 . However, if w = h u then the equilibrium point becomes unstable and a limit cycle with period T(τ * ) = 2π u appears there, see [35]. It follows from the numerical experiments that for τ > τ * , the Eq. (1) has only one stable cycle, and its period T(τ ) exceeds 2τ . So, we assume that this equation does not have other limit cycles with periods greater than 2τ , though we have not proved this analytically.
At the same time, this Eq. (1) can possess cycles whose periods are less than τ . Actually, for τ ≥ τ * and any natural m, there exists at least one limit cycle such that τ * + mT(τ * ) ≤ τ . This limit cycle is constructed as follows: let S m (τ ) = τ + mT(τ ), where T(τ ) is the period corresponding to the parameters α, h, τ . Numerical analysis of these cycles shows that this function grows monotonically with τ . When τ grows from τ * to τ , we have S m (τ * ) ≤ τ and S m (τ ) > τ . So, we see that if for some natural m the inequality τ * + mT(τ * ) ≤ τ holds, then the Eq. (1) with the parameters α, h, τ m has a limit cycle corresponding to the parameters α, h, τ .

Limit cycles of gene networks models for n > 1
Natural generalization of the model (1) to the case of circular gene network with n genetic elements ( Fig. 1e) has the form: are concentrations as above. For n = 2, the system of equations of the type (2) is the simplest model of artificial molecular trigger [36,37], for n = 3, this is a model of Elowitz-Leibler repressilator [38], see also [39,40] for more biological interpretations. Detailed mathematical analysis of delayed arguments phenomena in similar models is given in [41]. For its piece-wise version, the sufficient and necessary conditions of existence of a limit cycle were obtained in [42]. For similar 5-dimensional system, sufficient condition of existence of a cycle was obtained in [43], where some results on its stability were presented as well.
In our numerical experiments, we have seen that if the system (2) with the parameters α i , β i , h i , τ i has a limit cycle C 1 with coordinate functions (x 1 (t), x 2 (t), . . . , x n (t)), then the same system with the parameters α i , β i , h i ,τ i such that i=n i=1 τ i = i=n i=1τ i , has a limit cycle C 2 (y 1 (t), y 2 (t), . . . y n (t)) whose coordinate functions y i (t) coincide with x i (t) for all i up to some shift of t. So, one can assume here that i=n i=1 τ i < T, where T is the period of the cycles C 1 , C 2 , and that τ 1 = τ 2 = . . . = τ n .

The case of infinite value of the parameter h
When h i → ∞ for all i, then the system (2) reduces to the following one: where L i are the step functions: Similar systems were studied in [41][42][43] in analysis of circular gene networks, i = 2, 3, . . . , n.
The graphs of periodic solutions x j (t) of the system (3) have a standard form, see the Figure 2.
We consider now the symmetric simplified version of the system (3) in the case α 1 = α 2 = . . . = α n = α, β 1 = β 2 = . . . = β n = 1, in order to describe all periodic trajectories of the system (3) which are symmetric with respect to the cyclic permutations of the variables In this case this system is reduced to just one equation with constant delay τ : Its solution has the form: We assume here that then the Eq. (4) has a limit cycle C 1 with the period T, and such τ is unique on the interval (0, T/2). In this case characteristics of the graph on the Fig. 2 have the form Numerical analysis of this cycle shows that it is stable. Now, fix τ > 0, and α > 1, then the Eq. (4) has a limit cycle C 2 with the period T(τ ) such that This is unique limit cycle such that T(τ ) > 2τ . Surely, if this condition T(τ ) > 2τ is satisfied, then the cycles C 1 and C 2 coincide.
, then the Eq. (4) has a countable set of limit cycles.
Proof Let k ≥ 1 be a natural number. Consider τ k < τ such that T(τ k ) calculated by the formulae (6), (7) satisfies the condition τ k + kT(τ k ) = τ . Then the equation has a unique positive solution T(τ k ).
So, we have constructed a countable set of limit cycles {C k } such that, given k, the relations (8) determine the cycle C k , and the lemma is proved. Now, we give a description of a symmetric limit cycle of n-dimensional system (4) with the paraemter α. Let T be the period of this cycle, and (x 1 (t), x 2 (t), . . . , x n (t) be its coordinate T-periodic functions. Since the cycle is symmetric, we have Here, i = 1, 2, . . . , n−1, and positive natural p does not exceed (n−1). Thus, each positive solution of the equation corresponds to a limit cycle of the system (4), and given τ > 0, T(τ n,p 1 ,k 1 ) = T(τ n,p 2 ,k 2 ) if and only if k 1 = k 2 , and p 1 = p 2 .
For even values of n, no stable cycles have been discovered in our studies, both numerical and analytical.
Consider an inverse problem of identifications of parameters of 3-dimensional dynamical system where L i are as in (3), and α i > 1, i = 1, 2, 3. Piecewise linear systems of the type (10) in different dimensions were studied in gene networks modeling earlier, see for example [47][48][49]. As it was shown in [50], the system (10) has a unique limit cycle C with the periodic coordinate functions x(t), y(t), z(t); let T be its period, and [ 0, T] be a segment on the t-axis subdivided by the points 0 < t 1 < t 5 < t 2 < t 3 < t 4 < T so that the graph of the function y(t) is depicted on the Fig. 2. Unusual position of t 5 is chosen specially, in order to refer to this Figure. So, y 0 = y(0) = y(T) is the minimal value of y(t), y 2 = y(t 2 ) is its maximal value, and y( The shapes of the graphs of the functions x(t), z(t) are similar to that of y(t), though their maxima and minima are different as well, as the lengths of the segments in the partition of [ 0, T]. So, it follows from the formulae (5) that and these relations imply the following proposition:

Proposition 1 If the moments t j of extrema of the coordinate functions x(t), y(t), z(t)
are known from observations of the limit cycle C, then these maximal and minimal values, and the parameters α 1 , α 2 , α 3 can be expressed explicitly in terms of t j from the Eq. (11).
For example, Analogous inverse problems of identifications of parameters of higher-dimensional dynamical systems of the type (4) can be formulated and solved in a similar way, see for example, [51,52].

The case of finite values of h
Lemma 2 If the symmetric system has a symmetric T-periodic trajectory with coordinate functions (x 1 (t), . . . , x n (t)), then where p is an integer such that 0 ≤ p < n, and for each j, j = 1, 2, . . . , n, the function x j (t) is a T-periodic solution of the Eq. (1) with the parameters α, h,τ = τ + p n T. The converse is also true.
Thus, each symmetric limit cycle of the system (12) can be constructed from the periodic solution of (1) with parameters α, h,τ = τ + p n T, and we can formulate the following algorithm of construction of these symmetric limit cycles: Algorithm 1.
It follows from Lemma 2 that in order to construct symmetric limit cycles of the system (12), one has to find T-periodic solutions of the equation where p = 1, 2, . . . , (n − 1). Then for any limit cycle x(t) of the Eq. (14), the system (12) has a limit cycle (x 1 (t), x 2 (t), . . . , x n (t)) such that For autonomous sytems, when τ = 0, this algorithm has the following illustrative presentation. Consider the graph of the function y(τ ) = T(τ ) τ , and fix any integer p, such that 0 < p < n. Then for T(τ ) > 2τ , the limit cycle does exist, if the line y = n p intersects the graph of y(τ ), see the Fig. 3.
then there exists only one such cycle. If , then there exist two such cycles, Consider as examples the limit cycles of n-dimensional systems, for α = 1000, h = 10, n = 2, 3, . . . , 11. Each symmetric limit cycle of n-dimensional system (12) with the parameters α, h, and τ , can be obtained from T(τ )-periodic trajectories of the Eq. (1) with the parameters α, h, and τ .
For n = 2 there are no symmetric limit cycles, For n = 3 there exists one symmetric limit cycle with the phase lag T/3, Fig. 3a. For n = 4 there exist two symmetric limit cycles with the phase lag T/4, Figure 2 b. For n = 5 there exist two symmetric limit cycles with the phase lag T/5 and one limit cycle with the phase lag 2T/5, Fig. 3c.
For n = 6 there is one symmetric limit cycle with the phase lag T/3, and two limit cycles with the phase lag T/6, Fig. 3d.
For n = 7 there is one limit cycle with the phase lag 3T/7, one limit cycle with the phase lag 2T/7, and two limit cycles with the phase lag T/7, Fig. 3e.
For n = 8 there is one symmetric limit cycle with the phase lag 3T/8, and two limit cycles with the phase lag T/4, Fig. 3f.
For n = 9 there is one symmetric limit cycle with the phase lag 4T/9, and two limit cycles with the phase lag 2T/9, Fig. 3g.
For n = 10 there is one limit cycle with the phase lag 2T/5, one limit cycle with the phase lag 3T/10, and two limit cycles with the phase lag T 5 , Fig. 3h. For n = 11 there is one limit cycle with the phase lag 5T/11, one limit cycle with the phase lag 4T/11, two limit cycles with the phase lag 3T/11, and two limit cycles with the phase lag 2T/11, Fig. 3i.
Using the Algorithm 1, similar examples can be described for higher-dimensional versions of the system (12).

On non-uniqueness of limit cycles in models of circular gene networks
Another algorithm of construction of symmetric limit cycles of the system (11) is based on factorization of the dimension n, see [53].
Proposition 2 Let C n be a symmetric limit cycle of n-dimensional system (12), and x j (t) = x j+1 (t − pT n /n). This cycle C n is irreducible if and only if p and n are coprime.
If p and n are not coprime, their greatest common divisor equals m, and n = m · s, p = m · q, then C n is generated by the limit cycle C m of m-dimensional system of the form (12) such that x r (t) = x r+1 t − q T m m , r = 1, 2, . . . , m.

Corollary 1 Each symmetric limit cycle of the system (3) is either irreducible, or is generated by an irreducible cycle.
The limit cycle C n,s is called reducible, and if the system (12) has a symmetric limit cycle which cannot be constructed by the Algorithm 2, we call it irreducible.
Numerical experiments show that limit cycle C n is stable if and only if the dimension n is odd, n = 2m + 1, and x j (t) = x j+1 t − m T n n . Here j = 1, 2, . . . , n, as above.

Discussion
This paper is devoted basically to symmetric limit cycles in models of gene networks regulated by negative feedback loops only. Two algorithms of construction of these cycles are described. In this paper we have shown that 1. (a) Natural gene networks, not necessary circular, have oscillating potential since usually they contain circular gene subnetworks.
(b) Initiation of transcription/translation processes in these networks is realized by sufficiently powerful promotor/SD-site: α >> 1.
(d) For low-dimensional gene networks, the period between the initiation moment and the appearance of an active molecule in oscillations of their synthesis should exceed some critical value: τ > τ * .
2. At the same time, for n = 2, the systems (4) describe so called molecular triggers, and do not have limit cycles. Their analogues in the natural gene networks are genetic subsystems, that act as switches from one mode of operation to another under external conditions and/or regulatory signals [36,37].
3. Such molecular triggers can realize in perspective memory elements in design and construction of artificial genetic programs. The main attractive property of these molecular triggers is their bistability, and absence of limit cycles makes these systems more stable. As an example, the Fig. 4 shows a model of gene network composed by molecular triggers which realizes addition of binary numbers.
We assume here, that activity of the first genetic element of this trigger and inhibited state of the second one means "1", and vice versa, "0" corresponds to activity of the second element, and inhibited state of the first one. Figure 4b shows the Fragment 1 which calculates the sum of least lower-order digits of two binary numbers. Figure 4c corresponds to Fragment 2 which realizes summation of higher-order digits. In order to construct a gene network which adds two n-digit binary numbers, one has to join one Fragment 1 and (n − 1) Fragments 2, as it is shown on the Fig. 4d for n = 8. We assume here that different genetic elements of this network encode different proteins, so each regulatory binding is realised only by the proteins which correspond to the arrows of this Figure. As an illustration, consider the logical sum 10010010+11001001=101011011, see Fig. 4e. Here, the triggers {T 1 } should be arranged according to the digits of the first summand, and the triggers {T 2 }, respectively, to the digits of the second one. This can be done by corresponding protein addition to the medium. Then, the regulatory binding bring this gene network into a state where the triggers {S} desribe the required sum. Here, activity of the last trigger in {S} corresponds to the highest-order digit of this sum.
Modeling of similar gene networks and corresponding experiments can be used in construction of artificial genetic systems, in particular, in design of specialized genetic computers.
4. Nonstationary dynamics in models of circular gene networks with negative feedback mechanisms appears for high level of the Hill's parameter h. In the case of eukaryotes and prokaryotes, the necessary value of this nonlinearity parameter is achieved by multimerization of the transcriptional factors and/or by presence of several binding sites in DNA.
5. Nonstationary dynamics in these circuits is provided by the phenomenon of regulatory signal delay. In eukaryotes, this regulatory delay is achieved by cascade reactions in signaling pathways, separation of mRNA and protein synthesis into different compartments and transport of substances from one compartment to another. For one eukaryotic system, somitogenesis of vertebrata, as it is well-known, its development is regulated by an auto-oscillating process, controlled by autorepressor HES7 [24,[54][55][56]. Decreasing of the delay of autorepression signal of the gene HES7 as a result of introns deletion for this gene, leads to distortions of oscillation period, and even to their disappearance, if all the introns are removed [57,58].
Analysis of stability of oscillations in this system shows that the cooperativity coefficient h for HES1 self-repression should be sufficiently large (i.e. h ≥ 4) and the duration of the repression loop is between 40 and 60 min [59].
Thus, it is natural to suppose that for prokaryotes, due to the coupling of transcription and translation processes, the time period between the beginning of gene expresion and appearence of the active form of the transcriptional factor is too short for formation of oscillations.
Analysis of dynamics of synthetic constructions based on the prokaryotic genes and regulatory circuits with both, positive and negative loops shows that the regulatory signal delay that appears due to the negative autoregulatory loop is the key element of the oscillating construction design [60]. Moreover, numerical experiments show that extension of this loop leads to stable oscillations even in absence of positive regulatory loops [39].
From a theoretical viewpoint, the necessary value of the time lag can be realized by moving of the regulatory factor translation frame far enough from the transcription start. This can be done by long nucleotide sequence insertion between the promoter and transcription factor gene encoding. So, transcription of such long segments will imply the delay of appearence of the end product. Variations of the lengths of these insertions allow to control the time lag of the process.
This approach can be used in construction of synthetic gene networks with target functions. Corresponding mathematical models allow to find out structural, functional and parametric characteristics which provide necessary properties of the gene networks and give possibilities to plan strategy of their genetic synthesis.

Conclusions
We have proposed two algorithms of construction and localization of limit cycles in these gene networks models. Existence and stability of some of these cycles have been shown analytically [45], and for this gene networks model we have solved an inverse problem of parameters identification. Several natural questions on behavior of trajectories of the system (12) remain still open. In particular, stability of cycles of this system is not proved; here, we have just some results of numerical experiments.
The modeling results demonstrate that non-stationary dynamics in the models of circular gene networks with negative feedback loops is achieved by a high degree of non-linearity of the mechanism of the autorepressor influence on its own expression, by the presence of regulatory signal delay, the value of which must exceed a certain critical value, and transcription/translation should be initiated from a sufficiently strong promoter/Shine-Dalgarno site. We believe that the identified patterns are key elements of the oscillating construction design.
Existence of asymmetric limit cycles in the symmetric systems of the types (3) and (12) is also unclear. We plan to consider these questions in our forthcoming studies.

Methods
The main aim of our investigations is nonlocal analysis of some gene networks models represented by dynamical systems of kinetic type, and combimatorial structure of their phase portraits. In particular, we study here conditions of existence, stability, and (non)uniqueness of cycles of the models as well as periods of these cycles, and their geometric location in phase portraits of the dynamical systems.

Description of a model of simplest autorepressilator
The equation corresponds to a subsystem containing just one gene that encodes protein with concentration x = x(t) acting as the regulator of its own expression by repression mechanism (autorepressor). The parameters of this model (1) of the simplest circular regulatory circuit are interpreted as follows: the constant α describes efficiency of the gene expression, h is a Hill coefficient, it determines degree of non-linearity of positive monotonically decreasing function f, which describes negative feedback mechanism of gene expression regulation. The parameter τ is the lag in the delay of the argument, the subtrahend in this equation and in higher-dimensional models corresponds to velocity of natural loss in the network (degradation, dissipation etc).

Methods of analysis of higher-dimensional models
Our analysis of phase portraits of the higher-dimensional gene networks models with negative feedback loops, delayed arguments phenomena, and studies of corresponding inverse problems were based on classical methods of the qualitative theory of ordinary differential equations, and on our recent results [44,45] concerning monotonicity of the Poincaré map (the so-called succession map). Geometric arguments, topological fixed point theorem, and analytical calculations allow to prove stability of cycles and their nonuniqueness for some low-dimensional gene networks models.
At the same time, these geometric considerations show that phase portraits of ndimensional dynamical systems of the type (10) contain invariant domains which can be decomposed to unions of 2 n parallelepipeds (blocks) such that solutions of such systems have very simple analytic form in each block (Fig. 2). So, the parameters identification inverse problems for the systems (10) can be reduced to combinatorial description of transitions of their trajectories from block to block. In the case n = 3, the solution is given in Proposition 1. Higher-dimensional versions of this inverse problem require more combinatorial efforts even in the case n = 4, see [49].

Calculation method
The methods of integrating the systems of differential equations with delayed arguments are classical as well, see [28,30]. In our numerical experiments with these systems, we used the semi-implicit difference scheme with variable step. Corresponding calculations and construction of the graphs on the Figs. 2 and 3 which illustrate our mathematical results were carried out on personal computers of the first author.