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

Abstract

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.

Introduction

Molecular-genetic 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 sub-networks 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 X-ray sensitivity A) and MarR (multiple antibiotic resistance regulator) [4, 5], and for transcriptional factors H-NS (histone-like nucleoid structuring protein), IHF (integration host factor), Fis (factor for inversion stimulation), Lrp (leucine-responsive 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 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–18], etc.

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

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 sub-networks 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

$$ \frac{d x(t)}{dt} = f(x(t-\tau)) - x(t); \quad f(x) = \frac{\alpha}{1+x^{h}}. $$

(1)

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, 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₀. 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:

$$ \frac{d x_{1}}{dt} = \frac{\alpha_{1}}{1+x_{n}^{h_{n}} (t - \tau_{n})} - \beta_{1} x_{1}, \, \frac{d x_{i}}{dt} = \frac{\alpha_{i}}{1+x_{i-1}^{h_{i-1}} (t - \tau_{i})} - \beta_{i} x_{i}, \,\,\, i=2,3,\ldots n. $$

(2)

As in (1), the functions $f_{i}(x) = \frac {\alpha _{i}}{1+x^{h_{i}}}$ describe negative feedback mechanisms in the gene network, x₁=x₁(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 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₁ with coordinate functions (x₁(t),x₂(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₂(y₁(t),y₂(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₁,C₂, and that τ₁=τ₂=…=τ_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:

$$ \frac{d x_{1} (t)}{dt} = L_{1} (x_{n} (t-\tau_{1}))- \beta_{1} x_{1} (t); \quad \frac{d x_{i} (t)}{dt} = L_{i} (x_{i-1} (t-\tau_{i}))- \beta_{i} x_{i} (t); $$

(3)

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 α₁=α₂=…=α_n=α,β₁=β₂=…=β_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₁→x₂→…→x_n→x₁. In this case this system is reduced to just one equation with constant delay τ:

$$ \frac{dx}{dt} = \alpha - x(t) \,\,\,\, \text{for} \,\,\,\, x(t-\tau) < 1; \quad \frac{dx}{dt} = - x(t) \,\,\,\, \text{for} \, \,\,\, x(t-\tau) > 1; $$

(4)

Its solution has the form:

$$ \begin{aligned} & x(t) = (x(t_{i}) - \alpha) e^{-(t_{i} - t)} + \alpha \,\,\,\, \text{if} \,\,\,\, x(t_{i}) - t = 1 \,\,\,\, \text{and} \,\,\,\, x(t-\tau) < 1; \\ & x(t) = x(t_{i}) e^{-(t_{i} - t)} \,\,\,\, \text{if} \,\,\,\, x(t_{i} - t) = 1 \,\,\,\, \text{and} \,\,\,\, x(t-\tau) > 1. \end{aligned} $$

(5)

We assume here that 0<t₁<t₂<…<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₁ 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

$$ \begin{aligned} & t_{1}= \ln \frac{\alpha - e^{-\tau}}{\alpha - 1}, \quad t_{2}= \ln \frac{\alpha e^{\tau} - 1}{\alpha - 1}, \quad t_{3}=T-\tau, \,\, \\ & x_{0} =x(0)=x(T)=e^{-\tau}, \quad x(t_{2})= \frac{e^{\tau} \alpha - \alpha+1}{e^{\tau}}. \end{aligned} $$

(6)

Fig. 2

Standard form of the graphs of the solutions x_i(t) of the system (3). The horizontal axis is t, the vertical lines t₁,t₃ correspond to x_i(t)=1

Numerical analysis of this cycle shows that it is stable.

Now, fix τ>0, and α>1, then the Eq. (4) has a limit cycle C₂ with the period T(τ) such that

$$ T(\tau)= \ln \frac{(\alpha e^{\tau} - \alpha +1)(\alpha e^{\tau} -1)} {\alpha-1}. $$

(7)

This is unique limit cycle such that T(τ)>2τ. Surely, if this condition T(τ)>2τ is satisfied, then the cycles C₁ and C₂ 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

$$ e^{T(\tau_{k})} = \frac{\left(\alpha e^{\tau-kT(\tau_{k})} - \alpha +1\right) \left(\alpha e^{\tau - kT(\tau_{k})} -1\right)}{\alpha-1}; \quad \tau_{k} = \tau - kT(\tau_{k}) > 0 $$

(8)

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₁(t),x₂(t),…,x_n(t) be its coordinate T-periodic functions. Since the cycle is symmetric, we have

$$x_{i} (t) = x_{i+1} \left(t+ \frac{p}{n}T\right), \qquad x_{n} (t) = x_{1} \left(t+ \frac{p}{n}T\right). $$

Here, i=1,2,…,n−1, and positive natural p does not exceed (n−1). Thus, each positive solution of the equation

$$e^{T(\tau_{n,p,k})} = \frac{\left(\alpha e^{\tau + (p/n - k) T(\tau_{n,p,k})} -\alpha +1\right) \left(\alpha e^{\tau + (p/n - k) T(\tau_{n,p,k})} -1 \right)}{\alpha - 1}, $$

$$ \tau_{n,p,k} = \tau + \left(\frac{p}{n} -k \right) T(\tau_{n,p,k}) > 0, \,\,\, p= 1,2,\ldots n-1, \,\, k=0,1, \ldots \left[ \frac{\tau}{T(\tau_{n,p,k})} + \frac{p}{n} \right] $$

(9)

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₁=k₂, and p₁=p₂.

The case absence of delay arguments

In the particular case of the system (3) with τ₁=τ₂=…=τ_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, 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 3-dimensional dynamical system

$$ \frac{dx(t)}{dt}=L_{1} (z)-x(t); \,\, \frac{dy(t)}{dt}=L_{2} (x)-y(t); \,\, \frac{dz(t)}{dt}=L_{3} (y)-z(t), $$

(10)

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 t-axis subdivided by the points 0<t₁<t₅<t₂<t₃<t₄<T so that the graph of the function y(t) is depicted on the Fig. 2. Unusual position of t₅ is chosen specially, in order to refer to this Figure. So, y₀=y(0)=y(T) is the minimal value of y(t),y₂=y(t₂) is its maximal value, and y(t₁)=y(t₃)=1. Respectively, z₁=z(t₁) is the maximal value of z(t),z₃=z(t₃) is its minimal value, x₅=x(t₅) is the minimal value of x(t), and x₄=x(t₄) is its maximal value. At the same time, z(t₄)=z(t₅)=x(0)=x(t₂)=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

$$ \begin{aligned} & z_{1} \cdot e^{-(t_{5}-t_{1})} =1; \,\,\, z_{1} \cdot e^{-(t_{3}-t_{1})} = z_{3}; \,\,\, (z_{3} - \alpha_{3})\cdot e^{-(t_{4}-t_{3})} + \alpha_{3} = 1; \\ & y_{2} \cdot e^{-(t_{3}-t_{2})} = 1; \,\,\,\, y_{2} \cdot e^{-(T-t_{2}) }=y_{0}; \,\,\,\, (y_{0} - \alpha_{2}) \cdot e^{-t_{1}} + \alpha_{2} = 1; \\ & x_{4} \cdot e^{-(T-t_{4}) }= 1; \,\, x_{4} \cdot e^{-(T-t_{4})- t_{5}}= x_{5}; \,\, (x_{5} - \alpha_{1}) \cdot e^{-(t_{2}-t_{5}) }+\alpha_{1} = 1, \end{aligned} $$

(11)

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 α₁,α₂,α₃ can be expressed explicitly in terms of t_j from the Eq. (11).

For example,

$$z_{1} = e^{(t_{5}-t_{1}) }>1; \,\, z_{3} =e^{(t_{5}-t_{3})} <1; \,\, \alpha_{3} \left(1-e^{-(t_{4}-t_{3})}\right)= 1-e^{-(t_{4}-t_{5})} ; \quad \text{etc.} $$

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

$$ \frac{d x_{1}}{dt} = \frac{\alpha}{1+x_{n}^{h} (t - \tau)} - x_{1}, \,\, \ldots \, \frac{d x_{i}}{dt} = \frac{\alpha}{1+x_{i-1}^{h} (t - \tau)} - x_{i}, \,\, i=2,3,\ldots n ; $$

(12)

has a symmetric T-periodic trajectory with coordinate functions (x₁(t),…,x_n(t)), then

$$ x_{1} \left(t-\left(\tau+ \frac{p}{n}T \right) \right) = x_{n} (t), \ldots, x_{i} \left(t-\left(\tau+ \frac{p}{n}T \right) \right) = x_{i-1} (t), \quad i = 2,3,\ldots, n, $$

(13)

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 $\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 T-periodic solutions of the equation

$$ \frac{d x(t)}{dt} = f(x(t- (\tau + \frac{p}{n}T)) - x(t); \quad f(x) = \frac{\alpha}{1+x^{h}}, $$

(14)

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₁(t),x₂(t),…,x_n(t)) such that

$$x_{i} (t) = x\left(t - i \left(\tau + \frac{p}{n}T \right) \right). $$

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.

Fig. 3

Geometric interpretation of the Algorithm 1 for τ=0. The blue curve is the graph of the function $T\left (\frac {p}{n} \tau _{1} \right)$. The horisontal axis is t, the vertical one shows the values of $T\left (\frac {p}{n} \tau _{1} \right)$ and $\frac {n}{p}$. On each graph p=1,2,…,n−1. a n=3; b n=4; c n=5; d n=6; e n=7; f n=8; g n=9; h n=10; i n=11

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 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 $\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 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].

Algorithm 2.

Let n=s·ℓ, and s>1,ℓ>1 be such a factorization. If s-dimensional system (12) has a symmetric limit cycle C_s(x₁(t),…,x_s(t)) such that x_i(t)=x_i+1(t−T_s/s),i=1,…,s, then the n-dimensional system of the type (12) has a symmetric limit cycle C_n,s(x₁(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 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} \left (t-q \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/SD-site: α>>1.

(c) Autorepressor’s action is nonlinear: h>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.

Fig. 4

Hypothetic summating gene network; a memory cell; b Fragment 1 calculating the sum of the least lower order digits; c Fragment 2 calculating the sums of higher order digits; d summating gene network for n=8; e the states of the triggers in the summating process

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₁} should be arranged according to the digits of the first summand, and the triggers {T₂}, 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–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].

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 self-oscillating 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 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

$$\frac{d x(t)}{dt} = f(x(t-\tau)) - x(t); \quad f(x) = \frac{\alpha}{1+x^{h}}; $$

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 non-uniqueness for some low-dimensional gene networks models.

At the same time, these geometric considerations show that phase portraits of n-dimensional dynamical systems of the type (10) contain invariant domains which can be decomposed to unions of 2ⁿ 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.