Chaos and Hyperchaos in a Model of Ribosome Autocatalytic Synthesis

Any vital activities of the cell are based on the ribosomes, which not only provide the basic machinery for the synthesis of all proteins necessary for cell functioning during growth and division, but for biogenesis itself. From this point of view, ribosomes are self-replicating and autocatalytic structures. In current work we present an elementary model in which the autocatalytic synthesis of ribosomal RNA and proteins, as well as enzymes ensuring their degradation are described with two monotonically increasing functions. For certain parameter values, the model, consisting of one differential equation with delayed argument, demonstrates both stationary and oscillatory dynamics of the ribosomal protein synthesis, which can be chaotic and hyperchaotic dependent on the value of the delayed argument. The biological interpretation of the modeling results and parameter estimation suggest the feasibility of chaotic dynamics in molecular genetic systems of eukaryotes, which depends only on the internal characteristics of functioning of the translation system.

Ribosomes are key structural and functional elements of any cell. They not only provide the basic machinery for the synthesis of all proteins necessary for cell functioning during growth and division, but for biogenesis itself. For the last several decades, various control aspects of ribosome synthesis, assembly and degradation have been studied intensively in both prokaryotes and eukaryotes (see review [1][2][3][4][5][6][7][8][9][10][11] ). However, the dynamic aspects of functioning of this "protein factory" within a single cell cycle are still not clear.
From a biological point of view, ribosomes are self-replicating structures because they provide and control all stages of their own biogenesis. By participating in the synthesis of RNA polymerases, ribosomes assure the transcription of ribosomal and messenger RNAs and the subsequent translation of all structural and regulatory proteins that control maturation, assembly, functioning, recycling, degradation and synthesis of ribosomes de novo. That is why, reproduction of ribosomes is an autocatalytic process which at the molecular genetic level represents a gene network in which synthesis and degradation of rRNA, as well as ribosomal proteins and their regulators, are controlled by feedback mechanisms.
We have previously studied the role of retardation in generating chaos in genetic networks and have demonstrated that retardation existence (in an explicit or implicit form) is a prerequisite for chaos generation 17 , although quantitative estimates of retardation values are not known.
From the retardation point of view, transcription/translation systems in eukaryotic and prokaryotic organisms are contrasting genetic systems, yet main stages of their course are identical. The reason for this is the fact that prokaryotic transcription/translation processes are conjugated, that is, proceed almost simultaneously, while in eukaryotes they are separated into different compartments: transcription occurs in the nucleus and translationin the cytoplasm (Fig. 1).
Even in a simplified form demonstrated in Fig. 1, the complexity and duration of the eukaryotic ribosome reproduction, accompanied by multiple transport processes, compared to prokaryotic ribosome reproduction, are clearly visible. In this sense, one can say that retardation is a natural phenomenon in eukaryotic gene networks. Perhaps that is why we usually observe chaotic dynamics in mathematical models of natural eukaryotic gene networks 14,18,20-26 . However, it would be incorrect to a priori deny the existence of the retardation phenomenon in the gene networks of prokaryotic organisms. Gene networks of high dimensionality, controlled through complex nonlinear interactions, for which the formation of a chaotic dynamics was theoretically shown 27 , are characteristic of prokaryotes, likewise of eukaryotes.
From this point of view, the fundamental cellular processes common to all living organisms, including the biogenesis of ribosomes, which have chaos generating factors in their structural and functional organization, but demonstrate different dynamic characteristics in different organisms, are of particular interest.
In this paper, we present a theoretical analysis of the functioning dynamics using a model of a simple genetic system of ribosome biogenesis, which consists of two subsystems that describe positive autocatalytic synthesis of ribosomes and ribosome-degrading enzymes depending on the value of delay in the absence of any regulatory interactions.
We demonstrate that the model of ribosome biogenesis, in which the processes of ribosome synthesis and degradation are described with two monotonically increasing functions display chaotic and hyperchaotic dynamics of ribosome synthesis for certain parameter values. Previously, the presence of chaotic dynamics has been demonstrated only for systems described with unimodal and monotonically decreasing functions 15,[17][18][19]28 .

Model
The model describes a simple gene network of ribosome autocatalytic synthesis, which in general is the same in eukaryotic and prokaryotic organisms. By participating in the synthesis of RNA polymerases, ribosomes assure the transcription of ribosomal and messenger RNAs and the subsequent translation of regulatory proteins that control maturation, assembly and functioning of ribosomes. In this sense, ribosomes synthesize themselves. However, because RNA polymerase II synthesizes all kinds of messenger RNAs required for the cell functioning, it as well synthesizes the mRNA of protein-degrading enzymes (proteases) which have either lost their functionality, or the need for their presence in a cell has disappeared due to the changes in external conditions. By synthesizing these proteins, ribosomes enhance the process of their own degradation. Thus, the model describes positive processes of ribosome self-synthesis and ribosomal protease synthesis, the functional features of which lead to the negative ribosome number regulation. The genetic circuit corresponding to this system is depicted in Figure SI1.2.
The process of ribosome self-synthesis and ribosomal protease synthesis is described in the model with one equation with two monotonically increasing functions and one retarded argument. In this case, the retarded argument reflects the time taken by the cell for synthesis of ribosomal and messenger RNAs, their processing, precursor formations of small and large ribosomal subunits, their transport to the cytoplasm (in eukaryotes), maturation of ribosomal subunits, assembly of active ribosomes and mRNA translation of corresponding proteins.
The equation describing the considered gene network of ribosome autocatalytic synthesis (Fig. 2) is given below: Where, p(t) -intracellular concentration of ribosomes at the current time, f(x) -a control function of ribosome synthesis rate, g(x) -a dynamic parameter of the autodegradation rate, α -the rate constant for synthesis, β 1 -the rate constant for the constitutive degradation, β 2 -the rate constant for the dynamic autodegradation, K I -the efficiency constant for the autoactivation process, K D -the efficiency constant for the autodegradation process, h -the Hill coefficient that determines the nonlinearity degree of the dynamic autodegradation.
Justification of the functions f(x) and g(x) is given in Supplementary section SI1. In the following section we demonstrate that for certain values of the retarded argument τ in the equation (1) the chaotic dynamics is implemented with the following fixed set of dimensionless parameters: And for a given set of parameters we carry out a detailed analysis of the dynamic properties of the model (1). When selecting the set of parameter values, we proceeded from the fact that h I = 1 corresponds to the simple model of synthesis initiation and h D = 5 approximately corresponds to the stoichiometry of the prokaryotic proteases 29 . The values of other parameters when converted to the dimensional units are in the micromolar range of protein concentration and are charachteristic of cells (see Supplementary section SI6).
In the Supplementary section SI2 we provide an extensive analysis of the dependence of the chaotic potential of a simple model of ribosome biogenesis (1) on the parameters h I and h D , describing nonlinear processes of ribosome synthesis and degradation.
In the Supplementary section SI4 we provide a more complex version of the model of ribosome biogenesis, taking into account the multimeric architecture of ribosomes and the protein degradation system, as well as the two stage mechanism of translation initiation; and conduct, due to a high combinatorics, an analysis of the chaotic potential of the model (equations SI1, SI11, SI12) for the 13 selective sets of parameter values (Supplementary section SI5) depending on the delay parameter τ.

Results
Stationary solutions of the model (1), (2). In a model of ribosome biogenesis (1), the form of functions f(x) and xg(x), which control the rate of ribosome synthesis and autodegradation, suggests a zero steady state in the model (1) when x = 0. It is easy to determine, that the steady state is unstable when β 1 > 0. From the diagrams of the functions f(x) and xg(x) presented in Fig. 2, it follows that for the parameter values (2) the model represented by the equation (1) has a non-zero steady state solution. We find that for given values of the parameters (2) f(x) = xg(x) and is located at a ≅ . x 1826 1 0 point. According to Èl'sgol'ts & Norkin 29 , we determine the dependence of the stability by the first approximation of the found non-zero stationary solution of equation (1) on the parameter τ, τ ≥ 0. where Z(t) and Z(t − τ) are small deviations from the stationary solution.
After linearization of the right hand side of the equation (1) we obtain the linear equation with retarded argument in relation to Z(t): According to Èl'sgol'ts & Norkin 30 , implementation of this inequality means that when τ 0 > 0 there are such values of retardation parameter τ, at which a stationary solution of equation (1) is either asymptotically stable (τ < τ 0 ), or unstable (τ > τ 0 ). The value τ 0 can be found from the formula: It was calculated that τ 0~8 1.774 for the parameter values (2). That is, when τ ≥ 81.774, the equation (1), (2) has no stable steady states. Numerical analysis has shown that in this case solutions of the equation (1), (2) have an oscillatory dynamics that for certain values of the parameter τ may be periodic, quasi-periodic, or is unstable with respect to the initial data. Figure 3a demonstrates the presence of a periodic solution of equation (1), (2) when 470 ≤ τ ≤ 471.06. However, even at τ ≈ 471.07 bifurcation of a quasiperiodic solution occures, which have been shown to exist in the interval 471.07 ≤ τ ≤ 475.14 ( Fig. 3b).

Periodic solutions of the model represented by equation (1), (2).
A further increase of τ leads to the conversion of quasi-periodic solutions into periodic solutions at first of multiplicity 5 (Fig. 3d (1), (2) becomes unstable with respect to the initial data ( Fig. 4a) and the Poincaré map acquires complex fractal structure (Fig. 4c,e).
The bifurcation diagram (Fig. 5a) constructed from the intersection points of equation (1), (2) solutions and the Poincaré plane p(t) = 2500, p(t − τ) ≥ 1800 well demonstrates the transition to chaos when parameter τ value increases from 470 to 490; as well as it shows that the transition is carried out through a period doubling cascade of Feigenbaum scenario 31,32 .
We should note the unusual nature of the bifurcation diagram in the value ranges τ > 496. 25 and τ > 510 in which we observe an abrupt increase in the area of chaotic dynamics, as reflected in Fig. 5a, with rapid borders expansion of the bifurcation diagram on the ordinate axis (arrows 1 and 2).
Quantification of the solution instability with respect to the initial data revealed the presence of only one positive Lyapunov exponent in the solutions of the equation (1), (2) (Fig. 5b) when τ ≤ 500, indicating that a chaos found in the investigated area of the equation (1), (2) solutions is not a hyperchaos. However, when τ = 900, a hyperchaos is still detected in the equation (1), (2) (Fig. 5c, d). This is evidenced by the presence of two positive senior Lyapunov exponents in the area of the equation (1), (2) solutions (Fig. 5e). Computer experiment demonstrates that the transition to hyperchaos occures from chaos.

Biological Interpretation of the Modeling Results
The results presented above and in the Supplementary sections SI2 and SI5, were obtained for the equations in which the time, constants and variables were measured in dimensionless units. However, in order to evaluate the biological significance of the parameters for which the chaotic dynamics is observed in the studied models of ribosome biogenesis, it is necessary to assess what these parameter values may correspond to, if we turn to the natural units of measurement. When switching to the dimensional units, only values of Hill coefficients h I and h D do not change in the model, which are dimensionless within the meaning. To convert the dimensionless parameters of the model to dimensional parameters we used the translation algorithm, which is described in the Supplementary section SI3.   Using this algorithm, we calculated the values of parameters K I and K D measured in mkM, parameter α measured in mkM/h, parameters β 1 , β 2 measured in h −1 and parameter τ measured in hours, for all sets of parameter values for which the models displayed chaotic dynamics (see Supplementary sections SI6).
When estimating the parameters, we proceeded from known experimental data on the half-life of ribosomal proteins, which varies from 2 to 150 hours or more in higher and lower organizms 33,34 ; and total concentration of ribosomes in the eukaryotic cell, which we roughly estimated to be 16 μ M based on the average parameters of the eukaryotic cell volume (~500 μ m 3 ), the number of ribosomes in the cell (~5 × 10 6 ) and the ribosome molecular weight (4.2 MDa) [http://bionumbers.hms.harvard.edu/].
The fact that the free fraction of ribosomes, which can be 15% 35 , is exposed to degradation in the cell was also considered.
The parameter values for which the models of ribosome biogenesis displayed chaotic dynamics (Fig. 5а, Fig.  SI2.1 and SI5.1) are summarized in the . Given the relatively high structural complexity of the degradation system, in both prokaryotes and eukaryotes 29,36 , this limitation, in our view, is not critical.
With regard to the values of parameters K I , K D , which are the efficiency constants for the activation processes of ribosome synthesis and degradation, in all calculations they vary from 0.2 to 100 μ М for K I and from 2 to 70 μ М for K D , that is, they are in the micromolar range of protein concentration characteristic of any cell.
The rate constant for synthesis, α, was estimated to be ~4-32 μ M/h. The rate constant for the synthesis of ribosomes in the yeast cell was experimentally 37 estimated to be ~7 μ M/h, that is, of the same order as the one calculated from the model.
The delay parameter τ, as follows from the Table SI6.1, varies from 1.2 to 7 hours. This assessment was carried out for the 2 hours half-life of ribosomal proteins, that is, for the lower limit of the parameter value. And this means that 1.2 hour is the minimum delay at which the chaotic dynamics is realized. If the 10 hours half-life of ribosomal protein is used for the calculation of dimensional parameters, which is acceptable for both prokaryotes and eukaryotes 33,34 , the minimum delay value τ increases to 6 hours.
It is evident that the resulting estimates quite fit within the physiological borders of the parameters characteristic of molecular-genetic systems, except for the τ parameter, at least for prokaryotes.

Discussion
In this paper, the chaotic potential of a simple genetic system of ribosome biogenesis regulation during the cell cycle was analyzed using a mathematical model. The model describes positive processes of the ribosome synthesis by themselves and the synthesis of ribosome-degrading proteins (green and red ovals in Fig. 1) with two monotonically increasing function. The analysis of functional dynamics of the equation (1), (2) has revealed the possibility of both stationary and oscillatory dynamics formations in the system of ribosome biogenesis, that for certain values of the parameters and retarded arguments was sensitive to the minor changes in the initial data, that is, it was chaotic. We stress that this system has no positive-negative feedback regulatory loops.
Nevertheless, we explain the possibility of the chaotic dynamics generation by the fact that in the system of ribosome biogenesis this circuit is present in a latent form, as the ribosome degradation process simulates negative regulation and the autocatalytic synthesis -positive regulation. An analogy can be found between this observation and the genetic systems controlled by two regulatory loops, positive and negative, for which the chaotic dynamics is well studied [16][17][18][19] . However, the system of ribosome biogenesis is not identical to the gene networks with positive-negative regulation of gene expression efficiency and cannot be described with the respective models [16][17][18][19] .
It is important to note that when the values of the parameters in the model of ribosome biogenesis are transferred to the intrinsic dimensions, they lie in the range of values typical for the living systems. In this case, as follows from the Table SI6.1, the minimum estimate of the delay parameter τ , at which the chaotic dynamics is observed, varies from 1.2 to 7 hours. Accordingly, from the model it follows that prokaryotes most likely do not experience the chaotic potential of the system, since the ~1.2-6 hours retardation of the signal is too long for prokaryotic systems.
For eukaryotes, the conclusion is not so unambiguous. In this regard, we consider the somite morphogenetic differentiation in vertebrates, the amount of stages in which depends on the type of organism and can reach 300 or more, for example, in snakes 38,39 . The rhythmically repeated process of somite formation in vertebrates is controlled by molecular oscillator, and the duration of each somite pair formation depends on the organism type, on the somite position relative to the longitudinal axis, and in some organisms it depends on the external conditions of their development and can range from several minutes to several hours even in organisms with a reasonably fast ontogenesis 39 . Therefore, the retardation of ~1.5-7 hours can be quite realizable for eukaryotic organisms.
It should be noted that the considered models of ribosome synthesis and degradation do not impose any strict structural restrictions on the biological system. Indeed, the value h I = 1 indicates that the minimum single step model of ribosome synthesis initiation is sufficient for the formation of chaotic dynamics. The value h D = 4 is suitable for the whole range of valuesh I . This complexity is not critical for the protein degradation system. It is known that in the E.coli cell, the main proteolytic enzyme, Lon protease, is a hexamer and the eukaryotic proteasome includes tens of proteins 29,36 .
And if we consider that the increase of the complexity of the model of ribosome biogenesis (see Supplementary section SI5) leads to a decrease of the h D value until it reaches 2, at which the model (SI1, SI11, SI12) demonstrates chaotic dynamics, then one can assume that Hill coefficients of both simple and expanded models of ribosome biogenesis are consistent with prokaryotes and eukaryotes.
In this article we demonstrate the transition to chaos via a cascade of period-doubling bifurcations, that is, via the Feigenbaum scenario 31,32 . This scenario is one of the most commonly observed scenarios of transition to chaos and is well described in such biological systems as the environmental population models [40][41][42][43][44] ; the biophysical models of the electrical activity dynamics of neurons and pancreatic β -cells [45][46][47][48] ; the models of intracellular oscillations of calcium concentrations 49 and NADH 12 concentrations; as well as in the molecular genetic systemsthe model of alternative splicing 18 . Hence, in fact, this scenario is observed at all levels of life organization -from molecular to population level.
It should be noted that in the system of ribosome autocatalytic synthesis both chaos and hyperchaos are observed. The peculiarity of the hyperchaotic attractor, in contrast to the chaotic attractor, is the existence of more than one positive Lyapunov exponent (Fig. 5e). This is a more complex dynamic behavior of the system, which has been recently confirmed for the activity of the mollusc Clione limacina sensory neurons. According to the authors 50 , the hyperchaotic dynamics may be a mechanism that reflects its food-search strategy. But the reason why hyperchaos is more suitable for this purpose remains unclear.
In the presented model of ribosome biogenesis the hyperchaotic dynamics arises from the chaotic dynamics when retarded argument increases. From a biological point of view, the feasibility of this phenomenon at the intracellular level is not obvious, as in the model of molecular-genetic system of ribosome biogenesis hyperchaos is observed outside the physiological values of the retarded argument. However, the very fact of its discovery suggests such a possibility in the molecular-genetic systems for which such values of retarded argument may well be physiological.
From a mathematical point of view, the received result is original, since it demonstrates the presence of a chaotic potential in the system described with two monotonically increasing functions with one delayed argument. Earlier, similar results were shown for the systems described with unimodal functions with delayed argument and for the systems described with monotonically decreasing functions with two delayed arguments 15,[17][18][19]28 ; in the first case a chaotic potential was found for large Hill coefficients (≥ 10) 28 , and in the second -for the case of significant difference between the values of delayed arguments [17][18][19]28 , that seems an unlikely case for molecular-genetic systems.
In conclusion, it should be noted that the present work is the first to theoretically demonstrate the possibility of chaotic dynamics generation in the system of ribosome biogenesis in the biologically relevant areas of functioning of the eukaryotic cell, which depends only on the internal features of the translation system functioning, in the absence of any regulatory impacts. And that means that in the evolutionary aspect the cell faced the problem of metabolic instability rather early, even at the origin of the translation system. Nevertheless, at the molecular-genetic level of the biological organization, as opposed to the organism and population levels [51][52][53][54][55][56]  The lack of experimental evidence of the chaos generation at the intracellular level in vivo may indicate that during evolution the cell has found a solution to this problem -the possibility of stabilizing the system in the presence of chaos-generating factors. However, we do not know yet how the cell has solved the problem.

Methods
Calculation method. Integration of the equation was performed by the steps method 30 using the original program written in Fortran. Calculations were performed on a computing complex of the Information and Computing Center, Novosibirsk State University (http://www.nusc.ru).
Determination of the type of the oscillatory dynamics. The type of the oscillatory dynamics can be visualized by the characteristic Poincaré map: circular trajectory will generate the Poincare map that consists of a finite number of points over which the trajectory will go in a certain order; the quasicyclic trajectory is displayed on the plane as a combination of closed curves; the Poincaré map of a strange attractor is an infinite set of unordered points. Description of the Poincaré map construction is given below.
The criterion of the chaotic state of the oscillatory dynamics was considered to be the sensitivity to initial data of the equation (1) solution, which is detected by calculating the difference between the two trajectories of one variable, calculated with two identical models that start at zero time point with the initial functions that vary by a "small" value. If, as we iterate over values, the difference becomes comparable with the fluctuation amplitude, the dynamics is concluded to be sensitive to initial data.
To quantify the instability of the solution with respect to the initial data, the Lyapunov exponents were used, which allow to explicitly identify the motion pattern in the system. The calculation method is given below.
With the correct choice of the plane position, the Poincaré map gives us a clear idea of the type of oscillatory motion. The map is constructed as follows: we choose a part of the plane, which intersects the trajectory of the equation (1) solution at a certain point х 0 . The first point of the trajectory intersection with the selected part of the plane, which follows after the point х 0 , we shall designate х 1 = P(х 0 ). Then we define the second point of trajectory intersection with the selected part of the plane, х 2 = P(х 1 ), and so on. The set of intersection points is named the Poincaré map.
The Lyapunov exponents characterize the degree of exponential growth (or decay) of disturbances near the attractor trajectory. A spectrum of Lyapunov exponents is equal in number to the dimensionality of the phase space. For systems of equations with retarded arguments, for which the phase space is infinite, the number of Lyapunov exponents is also infinite. However, only a few senior Lyapunov exponents have a significant impact. The presence in the spectrum of at least one positive Lyapunov exponent is a measure of chaos. If a chaotic dynamics has two or more positive Lyapunov exponents, such a dynamics is identified to be hyperchaotic 57 .
Let us give a description of the method of calculation of several senior Lyapunov exponents for the retarded equation based on a modification of a known method of Benettin 58,59 . Let us consider the retarded equation: Let the initial distribution x(t) be given in the interval (− τ, 0), with a uniform partition with the step h N 0 Simultaneously with the equation (7) we consider the k number of perturbed equations are calculated, which include the norms of the vectors after orthogonalization, but before renormalization.
An evaluation of the Lyapunov exponents is obtained as follows: After a sufficiently large number of steps M the values λ k will converge. In our calculations, the mean values and standard deviations when Mτ > 10 6 were used to assess the senior Lyapunov exponents. We shall note that the solutions of retarded equations (7) are invariant with respect to the reference time, which corresponds to the disturbance which does not increase and does not decay, therefore, there must be at least one zero parameter in the spectrum of Lyapunov exponents.