 Research
 Open Access
 Published:
Limit cycles in models of circular gene networks regulated by negative feedback loops
BMC Bioinformatics volume 21, Article number: 255 (2020)
Abstract
Background
The regulatory feedback loops that present in structural and functional organization of moleculargenetic 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 nonstationary dynamics in the models of circular gene networks with negative feedback loops is achieved by a high degree of nonlinearity 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/ShineDalgarno site. We believe that the identified patterns are key elements of the oscillating construction design.
Introduction
Moleculargenetic systems or gene networks are molecular machines which decode genetic programs and bring them up to all the levels of living systems organization, from the molecular level (metabolic networks, signaling pathways, splicing, transport etc.) to subcellular, cellular, tissue, and organismic one. During the evolution, these systems have obtained ability to adequate reaction on the external signals and on changing of the environmental conditions. Their structure is based on regulatory circuits, i.e., on subnetworks composed by genes, which encode transcriptional factors and some types of enzymes (for example kinase), and their different combinations configure gene networks regulated by positive, negative, or variable feedback mechanisms.
From the evolutionary viewpoint, configuration of loops with negative feedback loops regulation is more simple task for a cell, than that of loops with positive feedback loops, because it can be accomplished by a simple blockage of the transcription initiation site by the regulatory protein without formation of complex structures that interact with the RNA polymerase.
This type of regulation is widespread in nature, and its simplest example is a subsystem containing just one gene which controls its own expression by negative feedback mechanism. This is the minimal regulatory circuit, it acts on itself (Fig. 1a,b), and there exists a great amount of such regulatory subsystems. For the prokaryotes, in particular for E.coli, this regulation mechanism is typical for global regulators, such as FNR (fumarate nitrate reduction regulator) and CRP (cAMP receptor protein) [1–3] for stress response regulators LexA (locus for Xray sensitivity A) and MarR (multiple antibiotic resistance regulator) [4, 5], and for transcriptional factors HNS (histonelike nucleoid structuring protein), IHF (integration host factor), Fis (factor for inversion stimulation), Lrp (leucineresponsive regulatory protein) etc [6–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 helixloophelix 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–18], etc.
The next stage of regulatory circuits’ complexity is represented by the twogene systems, where each gene encodes expression regulator of another gene (Fig. 1c). Here, the examples of negative regulation circuit are pairs of regulatory factors Sox3Snail and DDX21Snail. 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 noncoding 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 AMPKmTORNRF2 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–27] but also the possibility of complex, chaotic dynamics formation [28–31], including the hyperchaotic one [32–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.
Results
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 x_{0}>0, where \(x_{0}(1+x_{0}^{h})= \alpha \). 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 \frac {x_{0}^{h}}{1+x_{0}^{h}} >1, u^{2}=w^{2} 1\), and \( \tau \geq \tau ^{*}=\frac {\pi  \arctan (u)}{u}\) then the equilibrium point becomes unstable and a limit cycle with period \(T(\tau ^{*})=\frac {2 \pi }{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:
As in (1), the functions \( f_{i}(x) = \frac {\alpha _{i}}{1+x^{h_{i}}} \) describe negative feedback mechanisms in the gene network, x_{1}=x_{1}(t),x_{i}=x_{i}(t) 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 ElowitzLeibler 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 piecewise version, the sufficient and necessary conditions of existence of a limit cycle were obtained in [42]. For similar 5dimensional 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 \( \alpha _{i}, \beta _{i}, h_{i}, \bar {\tau }_{i} \) such that \( {\sum \nolimits }_{i=1}^{i=n} \tau _{i} = {\sum \nolimits }_{i=1}^{i=n} \bar {\tau }_{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 \( {\sum \nolimits }_{i=1}^{i=n} \tau _{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: L_{i}(x)=β_{i}α_{i} for 0≤x≤1; L_{i}(x)=0 for x>1. Similar systems were studied in [41–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 x_{1}→x_{2}→…→x_{n}→x_{1}. In this case this system is reduced to just one equation with constant delay τ:
Its solution has the form:
We assume here that 0<t_{1}<t_{2}<…<t_{k}<T is partition of the period [0,T] to intervals (t_{i},t_{i+1}) where either x(t−τ)<1, or x(t−τ)>1.
If \( \tau = \ln \frac {\alpha + \sqrt {\alpha ^{2} + 4(e^{T}  1)(\alpha 1)}} {2 \alpha }\) 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.
Lemma 1
If T(τ)<2τ, 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 ndimensional 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 Tperiodic 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 \( \tau > 0, T(\tau _{n,p_{1},k_{1}}) = T(\tau _{n,p_{2},k_{2}}) \) if and only if k_{1}=k_{2}, and p_{1}=p_{2}.
The case absence of delay arguments
In the particular case of the system (3) with τ_{1}=τ_{2}=…=τ_{n}=0, it was shown in [44], that for odd values of n such a system has a limit cycle if and only if α_{i}>1 for all i=1,2,…,n. For n=3 this cycle is unique and stable [45].
If n=5,β_{1}=β_{2}=…=β_{5}=1, all the functions L_{i} coincide, and \(2 \alpha >5+\sqrt {5}\), then the system (4) has at least two limit cycles, see [46].
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 3dimensional 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–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 taxis 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(t_{1})=y(t_{3})=1. Respectively, z_{1}=z(t_{1}) is the maximal value of z(t),z_{3}=z(t_{3}) is its minimal value, x_{5}=x(t_{5}) is the minimal value of x(t), and x_{4}=x(t_{4}) is its maximal value. At the same time, z(t_{4})=z(t_{5})=x(0)=x(t_{2})=1. 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 higherdimensional 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 Tperiodic 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 Tperiodic solution of the Eq. (1) with the parameters \( \alpha, h, \bar {\tau }=\tau + \frac {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 \( \alpha, h, \bar {\tau }=\tau + \frac {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 Tperiodic 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(\tau) = \frac {T(\tau)}{\tau }\), and fix any integer p, such that 0<p<n. Then for T(τ)>2τ, the limit cycle does exist, if the line \( y = \frac {n}{p}\) intersects the graph of y(τ), see the Fig. 3.
So, if \(\frac {n}{p} \leq 2\), or \( \frac {n}{p} > \max _{\tau > \tau ^{*}} \frac {T(\tau)}{\tau }\) then the system (12) does not have limit cycles whose coordinate funcitions satisfy the relations (13).
If \(\, \frac {n}{p} = \max _{\tau > \tau ^{*}} \frac {T(\tau)}{\tau }\), or \( 2 < \frac {n}{p} < \frac {T(\tau ^{*})}{\tau ^{*}} \), then there exists only one such cycle.
If \(\, \frac {T(\tau ^{*})}{\tau ^{*}} \leq \frac {n}{p} < \max _{\tau > \tau ^{*}} \frac {T(\tau)}{\tau }\), then there exist two such cycles,
Consider as examples the limit cycles of ndimensional systems, for α=1000,h=10,n=2,3,…,11. Each symmetric limit cycle of ndimensional 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 \( \frac {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 higherdimensional versions of the system (12).
On nonuniqueness 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].
Algorithm 2.
Let n=s·ℓ, and s>1,ℓ>1 be such a factorization. If sdimensional system (12) has a symmetric limit cycle C_{s}(x_{1}(t),…,x_{s}(t)) such that x_{i}(t)=x_{i+1}(t−T_{s}/s),i=1,…,s, then the ndimensional system of the type (12) has a symmetric limit cycle C_{n,s}(x_{1}(t),…x_{n}(t)) such that x_{j}(t)=x_{j+1}(t−T_{s}/s). We assume here, as usual, that if i=s, then i+1=1, respectively, if j=n, then j+1=1.
Proposition 2
Let C_{n} be a symmetric limit cycle of ndimensional 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 mdimensional system of the form (12) such that \(x_{r} (t)=x_{r+1} \left (tq \frac {T_{m}}{m} \right), r=1,2,\ldots,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} \left (t m \frac {T_{n}}{n} \right)\). 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/SDsite: α>>1.
(c) Autorepressor’s action is nonlinear: h>1.
(d) For lowdimensional 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 lowerorder digits of two binary numbers. Figure 4c corresponds to Fragment 2 which realizes summation of higherorder digits. In order to construct a gene network which adds two ndigit 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 highestorder 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 wellknown, its development is regulated by an autooscillating process, controlled by autorepressor HES7 [24, 54–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 selfrepression should be sufficiently large (i.e. h≥4) and the duration of the repression loop is between 40 and 60 min [59].
In the case of the prokaryotes, there is a great deal of examples of the negative autoregulation of gene expression efficiency [1–9], but no examples of stable selfoscillating functioning regimes in these autorepressilators are known.
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 nonstationary dynamics in the models of circular gene networks with negative feedback loops is achieved by a high degree of nonlinearity 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/ShineDalgarno 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 nonlinearity 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 higherdimensional models corresponds to velocity of natural loss in the network (degradation, dissipation etc).
Methods of analysis of higherdimensional models
Our analysis of phase portraits of the higherdimensional 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 socalled succession map). Geometric arguments, topological fixed point theorem, and analytical calculations allow to prove stability of cycles and their nonuniqueness for some lowdimensional 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. Higherdimensional 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 semiimplicit 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.
Availability of data and materials
All data generated or analyzed during this study are included in this published article.
Abbreviations
 AMPKmTORNRF2:

AMPK (adenosine monophosphateactivated protein kinase)
 mTOR:

Mammalian target of rapamycin
 NRF2:

Nuclear factor erythroid 2 related factor 2
 bHLH:

Basic helixloophelix factor
 cAMP:

Cyclic adenosine monophosphate
 ComK:

Competence transcription factor
 CRP:

Cyclic receptor protein
 Fis:

Factor for inversion stimulation
 FNR:

Fumarate nitrate reduction regulator
 HNS:

Histonelike nucleoid structuring protein
 HES:

Hairy and enhancer of split
 IHF:

Integration host factor
 LexA:

Locus for Xray sensitivity A
 Lrp:

Leucineresponsive regulatory protein
 MarR:

Multiple antibiotic resistance regulator
 RB:

Retinoblastoma
References
 1
Aiba H. Autoregulation of the Escherichia coli crp gene: CRP is a transcriptional repressor for its own gene. Cell. 1983; 32:141–9.
 2
Takahashi K, Hattori T, Nakanishi T, Nohno T, Fujita N, Ishihama A, Taniguchi S. Repression of in vitro transcription of the Escherichia coli fnr and nar X genes by FNR protein. FEBS Lett. 1994; 340:59–64.
 3
GonzálezGil G, Kahmann R, Muskhelishvili G. Regulation of crp transcription by oscillation between distinct nucleoprotein complexes. EMBO J. 1998; 17:2877–85.
 4
Brent R, Ptashne M. The lexA gene product represses its own promoter. Proc Natl Acad Sci USA. 1980; 77:1932–6.
 5
Prajapat MK, Jain K, Saini S. Control of MarRAB operon in Escherichia coli via autoactivation and autorepression. Biophys J. 2015; 109:1497–508.
 6
Falconi M, Higgins NP, Spurio R, Pon CL, Gualerzi CO. Expression of the gene encoding the major bacterial nucleotide protein HNS is subject to transcriptional autorepression. Mol Microbiol. 1993; 10:273–82.
 7
Aviv M, Giladi H, Schreiber G, Oppenheim AB, Glaser G. Expression of the genes coding for the Escherichia coli integration host factor are controlled by growth phase, rpoS, ppGpp and by autoregulation. Mol Microbiol. 1994; 14:1021–31.
 8
Wang Q, Wu J, Friedberg D, Plakto J, Calvo JM. Regulation of the Escherichia coli lrp gene. J Bacteriol. 1994; 176:1831–9.
 9
Pratt TS, Steiner T, Feldman LS, Walker KA, Osuna R. Deletion analysis of the fis promoter region in Escherichia coli: antagonistic effects of integration host factor and Fis. J Bacteriol. 1997; 179:6367–77.
 10
Shan B, Chang CY, Jones D, Lee WH. The transcription factor E2F1 mediates the autoregulation of RB gene expression. Mol Cell Biol. 1994; V(14):299–309.
 11
Hirata H, Yoshiura S, Ohtsuka T, Bessho Y, Harada T, Yoshikawa K, Kageyama R. Oscillatory expression of the bHLH factor Hes1 regulated by a negative feedback loop. Science. 2002; V(298):840–3.
 12
Bonev B, Stanley P, Papalopulu N. MicroRNA9 Modulates Hes1 ultradian oscillations by forming a doublenegative feedback loop. Cell Rep. 2012; 2:10–8.
 13
Navarro P, Festuccia N, Colby D, Gagliardi A, Mullin NP, Zhang W, KarwackiNeisius V, Osorno R, Kelly D, Robertson M, Chambers I. OCT4/SOX2independent Nanog autorepression modulates heterogeneous Nanog gene expression in mouse ES cells. EMBO J. 2012; 31:4547–62.
 14
Fidalgo M, Faiola F, Pereira CF, Ding J, Saunders A, Gingold J, Schaniel C, Lemischka IR, Silva JC, Wang J. Zfp281 mediates Nanog autorepression through recruitment of the NuRD complex and inhibits somatic cell reprogramming. Proc Natl Acad Sci USA. 2012; 109:16202–7.
 15
Trieu M, Ma A, Eng SR, Fedtsova N, Turner EE. Direct autoregulation and gene dosage compensation by POUdomain transcription factor Brn3a. Development. 2003; 130:111–21.
 16
Magenheim J, Hertz R, Berman I, Nousbeck J, BarTana J. Negative autoregulation of HNF4alpha gene expression by HNF4alpha1. Biochem J. 2005; 388:325–32.
 17
Monteiro R, Pouget C, Patient R. The gata1/pu.1 lineage fate paradigm varies between blood populations and is modulated by tif1 γ. EMBO J. 2011; 30:1093–103.
 18
Foka P, Singh NN, Salter RC, Ramji DP. The tumour necrosis factoralphamediated suppression of the CCAAT/enhancer binding proteinalpha gene transcription in hepatocytes involves inhibition of autoregulation. Int J Biochem Cell Biol. 2009; 41:1189–97.
 19
Zhang H, Zhang Y, Chen C, Zhu X, Zhang C, Xia Y, Zhao Y, Andrisani OM, Kong L. A doublenegative feedback loop between DEADbox protein DDX21 and Snail regulates epithelialmesenchymal transition and metastasis in breast cancer. Cancer Lett. 2018; 437:67–78.
 20
Gamba P, Jonker MJ, Hamoen LW. A novel feedback loop that controls bimodal expression of genetic competence. PLoS Genet. 2015; 11:e1005047.
 21
Delihas N. Regulating the regulator: MicF RNA controls expression of the global regulator Lrp. Mol Microbiol. 2012; 84:401–4.
 22
Holmqvist E, Unoson C, Reimeg ard J, Wagner EG. A mixed double negative feedback loop between the sRNA MicF and the global regulator Lrp. Mol Microbiol. 2012; 84:414–27.
 23
Kapuy O, Papp D, Vellai T, Bánhegyi G, Korcsmáros T. Systemslevel feedbacks of NRF2 controlling autophagy upon oxidative stress response. Antioxidants (Basel). 2018; 7:39.
 24
Kageyama R, Niwa Y, Isomura A. González A, Harima Y. Oscillatory gene expression and somitogenesis. Wiley Interdiscip Rev Dev Biol. 2012; 1:629–41.
 25
Tsai TY, Choi YS, Ma W, Pomerening JR, Tang C, Ferrell JR. Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science. 2008; 321:126–9.
 26
Shen S, Ma Y, Ren Y, Wei D. Construction of an oscillator gene circuit by negative and positive feedbacks. J Microbiol Biotechnol. 2016; 26:139–44.
 27
Edwards R, Glass L. Combinatorial explosion in model gene networks. Chaos. 2000; 10:691–704.
 28
Khlebodarova TM, Kogai VV, Fadeev SI, Likhoshvai VA. Chaos and hyperchaos in simple gene network with negative feedback and time delays. J Bioinform Comput Biol. 2017; 15:1650042.
 29
Suzuki Y, Lu M, BenJacob E, Onuchic JN. Periodic, quasiperiodic and chaotic dynamics in simple gene elements with time delays. Sci Rep. 2016; 6:21037.
 30
Kogai VV, Likhoshvai VA, Fadeev SI, Khlebodarova TM. Multiple scenarios of transition to chaos in the alternative splicing model. Int J Bifurcat Chaos. 2017; 27:1730006.
 31
Edwards R, Siegelmann HT, Aziza K, Glass L. Symbolic dynamics and computation in model gene networks. Chaos. 2001; 11:160–9.
 32
Likhoshvai VA, Kogai VV, Fadeev SI, Khlebodarova TM. Chaos and hyperchaos in a model of ribosome autocatalytic synthesis. Sci Rep. 2016; 6:38870.
 33
Likhoshvai VA, Fadeev SI, Kogai VV, Khlebodarova TM. On the chaos in gene networks. J Bioinform Comput Biol. 2013; 11:1340009.
 34
Likhoshvai VA, Kogai VV, Fadeev SI, Khlebodarova TM. Alternative splicing can lead to chaos. J Bioinform Comput Biol. 2015; 13:1540003.
 35
Fadeev SI, Likhoshvai VA. On hypothetical gene networks. J Appl Indust Math Sib Zh Ind Mat. 2003; 6:134–53. (Russian).
 36
Gardner TS, Cantor CR, Collins JJ. Construction of a genetic toggle switch in Escherichia coli. Nature. 2000; 403:339–42.
 37
Tchuraev RN, Stupak IV, Tropynina TS, Stupak EE. Epigenes: design and construction of new hereditary units. FEBS Lett. 2000; 486:200–2.
 38
Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature. 2000; 403:335–8.
 39
Maeda K, Kurata H. Long negative feedback loop enhances period tunability of biological oscillators. J Theor Biol. 2018; 440:21–31. https://doi.org/10.1016/j.jtbi.2017.12.014.
 40
Sun M, Cheng X, Socolar JE. Causal structure of oscillations in gene regulatory networks: Boolean analysis of ordinary differential equation attractors. Chaos. 2013; 23(2):025104. https://doi.org/10.1063/1.4807733.
 41
Glyzin SD, Kolesov AY, Rozov NK. QuasiStable Structures in Circular Gene Networks. Comput Math Math Phys. 2018; 58(5):659–79.
 42
Golubyatnikov VP, Minushkina LS. Monotonicity of the Poincaré mapping in some models of circular gene networks. J Appl Ind Math. 2019; 13:472. https://doi.org/10.1134/S1990478919030086.
 43
Golubyatnikov VP, Golubyatnikov IV, Likhoshvai VA. On the Existence and stability of cycles in fivedimensional Models of Gene Networks. Numer Anal Appl. 2010; 3:329–35.
 44
Golubyatnikov VP, Ivanov VV. Cycles in odddimensional models of circular gene networks. J Appl Ind Mat. 2018; 12:648–57.
 45
Golubyatnikov VP, Ivanov VV. Uniqueness and stability of a cycle in 3dimensional blocklinear circular gene network models. Sib J Pure Appl Math. 2018; 18:19–28. (Russian).
 46
Akinshin AA, Golubyatnikov VP, Golubyatnikov IV. On some multidimensional models of gene network functioning. J Appl Ind Mat. 2013; 7:1–7.
 47
Glass L, Pasternack JS. Stable oscillations in mathematcal models of biological control systems. J Math Biol. 1978; 6:207–23.
 48
De Jong H, Gouze JL, Hernandez C, Page M, Sari T, Geiselmann J. Qualitative simulation of genetic regulatory networks using piecewiselinear models. Bull Math Biol. 2004; 66:301–40.
 49
Ayupova NB, Golubyatnikov VP. On two classes of nonlinear dynamical systems: The fourdimensional case. Sib Mat J. 2015; 56:231–6.
 50
Ayupova NB, Golubyatnikov VP. On the uniqueness of a cycle in an asymmetric threedimensional model of a melecular repressilator. J Appl Ind Mat. 2014; 8:1–6.
 51
Tieri P, Farina L, Petti M, Astolfi L, Paci P, Castiglione F. Network Inference and Reconstruction in Bioinformatics. Reference Module in Life Sciences. Encycl Bioinforma Comput Biol. 2019; 2:805–813. https://doi.org/10.1016/B9780128096338.202902.
 52
Delgado FM, GómezVela F. Computational methods for Gene Regulatory Networks reconstruction and analysis: A review. Artif Intell Med. 2019; 95:133–45. https://doi.org/10.1016/j.artmed.2018.10.006.
 53
Akinshin AA, Golubyatnikov VP. On cycles in symmetric dynamical systems. Bull Novosibirsk State Univ. 2012; 12:3–12. (Russian).
 54
Hirata H, Bessho Y, Kokubu H, Masamizu Y, Yamada S, Lewis J, Kageyama R. Instability of Hes7 protein is crucial for the somite segmentation clock. Nat Genet. 2004; 36:750–4.
 55
Mara A, Holley SA. Oscillators and the emergence of tissue organization during zebrafish somitogenesis. Trend Cell Biol. 2007; 17:593–9.
 56
Harima Y, Kageyama R. Oscillatory links of Fgf signaling and Hes7 in the segmentation clock. Curr Opin Genet Dev. 2013; 23:484–90.
 57
Takashima Y, Ohtsuka T, Gonzalez A, Miyachi H, Kageyama R. Intronic delay is essential for oscillatory expression in the segmentation clock. Proc Natl Acad Sci USA. 2011; 108:3300–5.
 58
Harima Y, Takashima Y, Ueda Y, Ohtsuka T, Kageyama R. Accelerating the tempo of the segmentation clock by reducing the number of introns in the Hes7 gene. Cell Rep. 2012; 3:1–7.
 59
Bernard S, Cajavec B, PujoMenjouet L, Mackey MC, Herzel H. Modelling transcriptional feedback loops: the role of Gro/TLE1 in Hes1 oscillations. Philos Trans A Math Phys Eng Sci. 2006; 364:1155–70.
 60
Stricker J, Cookson S, Bennett MR, Mather WH, Tsimring LS, Hasty J. A fast, robust and tunable synthetic gene oscillator. Nature. 2008; 456:516–9.
Acknowledgements
Not applicable.
About this supplement
This article has been published as part of BMC Bioinformatics Volume 21 Supplement 11, 2020: Selected Topics in “Systems Biology and Bioinformatics” 2019: bioinformatics. The full contents of the supplement are available online at  https://bmcbioinformatics.biomedcentral.com/articles/supplements/volume21supplement11.
Funding
Publication costs have been funded by the Program for Fundamental Research of the Siberian Branch RAS (project 032420190040C01, VAL and TMK) in the design of the study and collection, in interpretation of data, and in writing the manuscript, and by Russian Foundation for Basic Research (project 180100057, VPG) in analysis, and in writing the manuscript.
Author information
Affiliations
Contributions
VAL and TMK conceived the research. VAL and VPG carried out the theoretical proofs. VAL, VPG and TMK interpreted the results and wrote the paper. All authors edited and approved the final version of the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interest.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Likhoshvai, V.A., Golubyatnikov, V.P. & Khlebodarova, T.M. Limit cycles in models of circular gene networks regulated by negative feedback loops. BMC Bioinformatics 21, 255 (2020). https://doi.org/10.1186/s1285902003598z
Received:
Accepted:
Published:
Keywords
 Mathematical modeling
 Circular gene networks
 Delay argument equations
 Feedback loops regulation
 Autorepressor
 Cycles
 Phase portraits
 Inverse problems