Self-oScillatory dynamicS of the metabolic proceSS in a cell

In this work, a mathematical model of self-oscillatory dynamics of the metabolism in a cell is studied. the full phase-parametric characteristics of variations of the form of attractors depending on the dissipation of a kinetic membrane potential are calculated. the bifurcations and the scenarios of the transitions «orderchaos», «chaos-order» and «order-order» are found. We constructed the projections of the multidimensional phase portraits of attractors, Poincaré sections, and Poincaré maps. the process of self-organization of regular attractors through the formation torus was investigated. the total spectra of Lyapunov exponents and the divergences characterizing a structural stability of the determined attractors are calculated. the results obtained demonstrate the possibility of the application of classical tools of nonlinear dynamics to the study of the selforganization and the appearance of a chaos in the metabolic process in a cells.

м а т е м а т и ч н е мод е л юв а н н я біох і м і ч н и х п р оц е с і в м а т е м а т и ч н е мод е л юв а н н я біох і м і ч н и х п р оц е с і в T he development of a united theory of selforganization of the matter is one of the most important problems of natural sciences. Its solution is possible only at the comprehension of the unity of the physical laws for alive and inorganic matters. It is worth noting the significant contribution of synergetics to this field of knowledge. The idea of the appearance of dissipative structures in any open nonlinear systems is a corner stone in the foundation of such a theory. Synergetics explains the formation of an order from a chaos and, by this, outlines a way to the description of the biochemical evolution of the inorganic matter to the alive one [1][2][3][4].
The development of a general mathematical model of cell is a very complicated task. To solve it, it is necessary to study the biochemical processes and the types of nonlinearities characteristic of enzyme-substrate interactions causing the appearance of the self-organization in a system. Namely this problem belongs to the field of interests of the author.
Here, we deal with a mathematical model of the metabolism running in a cell arthrobacter globiformis. This model of enzyme kinetics is con-structed on the basis of the original method developed by V. P. Gachok [5,6] who proposed to apply a function of the general form V(X ) = aX/(1 + bX ) describing the adsorption of an enzyme in the domain of local coupling with a substrate to the modeling of biochemical processes.
As was noticed by a number of researchers [7,8], Prof. V. P. Gachok unified the construction of models for a wide class of biochemical systems described with multidimensional nonlinear differential equations.
In the present work, we will consider the meta bolic process in the arthrobacter globiformis cells quite well studied in experiments [9]. Bacteria of the given species are widely spread in the Nature and intensively used as producers of various metabolites. These bacteria are applied to the purification of wastes of the petrochemical, cokechemical, tannic, etc. productions, because they are able to decompose practically all hydrocarbons contained in oil. They are used to eliminate the toxic action of herbicides on plants and in other fields of ecology. In addition, these bacteria are involved in medicine in the diagnostics of the galactose exchange and in the transformation of steroids. Thus, the significance of a. globiformis culture for the industry is obvious. It is also worth noting that its cultivation in the continuous mode is characterized by the appearance of oscillatory phenomena in the metabolism of the given cells [10,11], which should be taken into account in their usage.
Such oscillatory mode was theoretically forecasted in work [12] prior to the experiment. Analogous modes of metabolic self-oscillations were registered in the processes of photosynthesis, glycolysis, variation of calcium in cells, oscillations in heart muscle, etc. [13][14][15][16]. Namely the oscillatory modes allow one to study various structural-functional couplings, which lead to the self-organization of metabolic processes in cells and a whole biosystem.
As distinct from the earlier results obtained within the model under consideration, we will calculate phase-parametric characteristics of variations of the form of attractors depending on the dissipation of a kinetic membrane potential and will show the appearance of the mutual transitions between stable and chaotic oscillations. With the help of a mathematical apparatus of the nonli near dynamics, we will perform the investigation of these transitions in the mentioned modes. The calculated scenarios such as «order-chaos», «chaosorder» and «order-order» and their mathematical analysis give a certain contribution to the understanding of the self-organization of biochemical systems.

a mathematical model and methods of its study
The general scheme of metabolic process running in arthrobacter globiformis cells under the transformation of steroids is presented in Fig. 1. The input substances for a bioreactor are hydrocortisone and oxygen. The output products are prednisolone and its hydroxy derivative. The process occurs in the following way.
Hydrocortisone (G) present in the flowing solution induces the biosynthesis of 3-ketosteroid-Δ′dehydrogenase enzyme (e 1 ) on the external side of a cytoplasmatic membrane of cells. Under the action of this enzyme, hydrocortisone is transformed into prednisolone (P). The 3-ketosteroid-Δ′dehydrogenase enzyme acquires the reducing form (e 1 ). By reducing the respiratory chain (q), this en-zyme is oxidized into (e 1 ). Oxygen (o 2 ) from the flowing solution oxidizes the respiratory chain into the form (Q). Produced prednisolone induces the biosynthesis of 20β-oxysteroid-dehydrogenase enzyme (e 2 ), which together with NaD⋅H(N ) coenzyme produces 20β-oxyderivative of prednisolone (B). Then a part of produced 20β-oxyderivative of prednisolone is consumed by cells in the Krebs cycle, by increasing the level of NaD⋅H. In the enhanced concentration, NaD⋅H blocks the respiratory chain. In this case, the process of transformation of hydrocortisone becomes slow.
The transformation of steroids is accompanied by the formation of a kinetic membrane potential (ψ). Its increased level holds the respiratory chain in the reduced state (q). The membrane potential is also used by cells in the Krebs cycle.
Here, V(X ) = X/(1 + X ); V (1) (ψ) = 1/(1+ψ 2 ); V(X) is a function related to the adsorption of the enzyme in the domain of local coupling (it was proposed by V. P. Gachok [5,6]); and V (1) (ψ) is a function characterizing the influence of the kinetic membrane potential on the respiratory chain.
This autonomous system of nonlinear differential equations was numerically solved by the Runge-Kutta-Merson method. The accuracy of calculations was set to be 10 -8 . To ensure the reliability of the study, namely the transition of the system from the initial and transient phases to the asymptotic state with an attractor, the computation time was taken to be 1 000 000. For this time, the trajectory was «stuck» onto the corresponding attractor.
The given system is dissipative. It is characterized by the dissipation of the concentrations of its reagents (except for the respiratory chain). Therefore, the elements of its phase space shrink during the evolution. Let us denote the set of all initial points of the system by V: , let as t → +∞ and L∈V. The limiting set L is an attracting one, and the system approaches an attractor, so that V is called the basin of a pulling attractor. We have obtained the following attractors: stationary modes (stable focus and stable node), self-oscillatory modes with various multiplicities and strange attractors.
We analyzed the stability of modes of the metabolic process with the use of the theory of stability of solutions of differential equations [31,32], which was already applied to the study of manystage enzymic reactions [33,34].
For the unambiguous identification of the type of an attractor, we calculated the total spectrum of Lyapunov exponents, for chosen points, λ 1 , λ 2 ,…,λ 10 and determined their sum: (see Table).
The calculation was carried out by Benettin's algorithm. The orthogonalization of the perturbation vectors was performed within the Gram-Schmidt method [31,35,36].
The studies within the model were carried out by the construction of phase-parametric characteristics, by changing one of the parameters (see below).
In order to identify the structure of attractors at the chosen points, we used the method of Poincaré sections. We note that the Poincaré section possesses properties of a flow and can be studied easier than the flow itself. In the phase space, we drew the cutting plane P = 0.2 and determined the coordinates of intersections for the corresponding variables on the two-dimensional projections of the constructed section. By the positions of the intersection points in this plane, we can judge about the structure of an attractor and the topological configuration of its phase trajectory.
Using the obtained Poincaré sections, we constructed the Poincaré maps, which allow one to determine the stationary mode of oscillations in the system.

results of studies
System (1) -(10) is a completely determined nonlinear system of differential equations in the ten-dimensional phase space. It describes the biochemical system with positive feedback, which generates the autocatalysis of NaD⋅H (see Eqs. (3) and (9)). Therefore, the self-oscillatory modes can appear in it. With the help of numerical computations, we have studied all possible modes arising in the given system at a variation of the dissipation of the kinetic membrane potential α (10). This potential is the most general parameter determining the stability of the metabolic process in a cell.
In Fig. 2,a, we present a phase-parametric characteristic of the system G(α) for the interval, where α∈(0.032,0.033). To construct the phaseparametric characteristic, we used the method of sections. In the phase space with a trajectory of the system, we drew the cutting plane P = 0.2. Such a choice is supported by the symmetry of oscillations relative to this point in multiple modes. As an example, Fig. 2,b,c shows the kinetics of G(t) and P(t) and the projections of phase portraits, which are used in the construction of the phase-parametric characteristic. When the trajectory P(t) intersects this plane in some single direc-tion, we mark the value of chosen variable (e.g., G) on the phase-parametric characteristic. In such a way, we determine the point of the intersection of the trajectory with the two-dimensional plane. If a multiple periodic limiting cycle arises, we obtain a number of points on the plane, which are repeated in the period. In the case where a deterministic chaos appears, the points of intersection of the trajectory with the plane are located chaotically.
In Fig. 2,a, we present the change of the multiplicity of oscillations in various modes for the whole interval of variation of α. There, the regions of the formation of regular and chaotic attractors are clearly seen. By examining the figure from right to left, i.e., by decreasing α from 0.033 down to 0.032, we see the appearance of the sequence of regular n⋅2 0 and strange n⋅2 x attractors, whose cycle period corresponds to the multiplicity n. If the given attractor is formed on a toroidal surface, it is denoted by the mark (t) (Table).
The less the coefficient of dissipation of the kinetic membrane potential of a cell, the more the energy is conserved in it, which complicates its internal dynamics. The 8-fold periodic cycle of a regular attractor 8⋅2 0 seen at α = 0.033 is broken with the formation of a strange attractor 8⋅2 x . In other words, the transition «order-chaos» occurs. The further decrease in α causes the self-organization of the system and the formation of a 9-fold periodic cycle, which corresponds to the transition «chaos-order». The given regular attractor 9⋅2 0 is total spectra of Lyapunov exponents for attractors of the system under study (λ 4 -λ 9 are not important for our investigation)  replaced by a strange attractor 9⋅2 x , then by 10⋅2 0 , 10⋅2 x , and so on up to the regular attractor of a quasiperiodic cycle ≈ n⋅2 0 . On this section of the phase-parametric characteristic, the decrease in the coefficient of dissipation of the kinetic membrane potential is periodically accompanied by the regular self-organization of the system. There occur the destruction of simple structures and the appearance of more complex ones. But, after the attractor ≈ n⋅2 0 , a simple regular attractor 1⋅2 0 appears unexpectedly. Let us consider the specific features of the self-organization of the system on separate intervals of the phase-parametric trajectory in more details. In Fig. 2,b,c, we present the kinetic behavior of the variables G(t) and P(t) and the projections of phase portraits for α = 0.032 and α = 0.0321107, respectively. In the interval α∈(0.032,0.03211), the system transits from a simple 1-fold periodic cycle to a 1-fold periodic cycle on a torus (Table). In this case, the amplitude increases.
At α = 0.0321107, a simple 8-fold periodic cycle appears instantly from the 1-fold cycle without intermediate doubling. Its kinetics and phase portrait are shown in Fig. 2,с. In other words, we observe the instantaneous transition «order-order». The found cycles are stable by Lyapunov, because their first Lyapunov exponents are equal to zero (Table).
By analyzing this transition, we note that the divergence of the system Λ is negative at the intermediate point α = 0.03211 and is larger in modulus than that at α = 0.032 and 0.0321107. For this mode, λ 1 , λ 2 , and λ 3 are zero. But, for the modes with α = 0.032 and 0.0321107, λ 2 and λ 3 are negative. This means that the element of the phase volume of the system at these points shrinks, as a whole, more strongly along the trajectory, i.e., the trajectory moves away from the toroidal surface (α = 0.03211) onto simple regular attractors.
The increase of |Λ| at α = 0.03211 is regular, because the dissipation (α) increases. As α increases up to 0.0321107, |Λ| decreases. This means that, after the appearance of a bifurcation and the transition from 1⋅2 0 to 8⋅2 0 , the metabolism in a cell organizes itself so that the total dissipation of the cell energy |Λ| becomes less at a higher dissipation of the kinetic membrane potential α. We note that the energy dissipation under oscillations realized on a toroidal surface in the multidimensional space is higher than that under simple periodic oscillations.
As α increases, the formed cycles are slowly shifted, on the whole, in the phase space in the single direction. It is seen in Fig. 2,a from inclined displacements of the traces of intersection points.
Let us now consider each part of the phaseparametric characteristic in more details. In Fig. 3,a, we show a part of the phaseparametric characteristic for α∈(0.03211,0.03215). As α increases to 0.03212, we observe the alternation of attractors. For example, a one-fold quasiperiodic cycle ≈ 1⋅2 0 arises again on a toroidal surface at α = 0.0321125. Then this cycle is broken at α = 0.03211295, and a 8-fold strange attractor 8⋅2 x arises. For this attractor, a part of the projection of its phase portrait is shown in Fig. 3,b. At α = 0.03212, the oscillations become simple, and an ordinary 14-fold periodic cycle is established (Fig. 3,с). As α increases to 0.0321286, the oscillations transit onto a toroidal surface. Quasiperiodic oscillations ≈ n⋅2 0 with irrational ratios of their frequencies are observed. The initial stage of the filling of the phase space by the trajectory is demon strated in Fig. 3,d by a part of the projection of the phase portrait.
As α increases further to 0.032135 (Fig. 3,a), we see the formation of lacunes on the phaseparametric characteristic, which means that the given multiple quasiperiodic attractors ≈ n⋅2 0 alternate with simple periodic 14-fold attractors. At α = 0.032135 as a result of the superposition of two periodic multiple oscillations with rationally commensurable frequencies, the strictly periodic mode 47⋅2 0 of oscillations on a torus is established. Its part is shown in Fig. 3,e.
In Fig. 4,a, we present the kinetic curve of this mode for the variable e 1 (t) (most clearly pronounced). The plot shows the periodic repetition of stable 47-fold oscillations.
The further increase in the parameter α leads to a change in the periodicity of oscillations on the toroidal surface. The quasiperiodic oscillations ≈39⋅2 0 are established at α = 0.032136, strictly periodic oscillations 39⋅2 0 at α = 0.032137, and then periodic 70-fold oscillations at α = 0.032138. At α = 0.032139, the oscillations become simple 14-fold. Their kinetics for the variable e 1 (t) (α = 0.0321395) is shown in Fig. 4,b, and the projection of the phase portrait of the regular attractor 14⋅2 0 is analogous to that in Fig. 3,с. On the considered interval of variation of α, we observe the alternation of attractors with different complexities and the transitions: «orderorder», «order-chaos», and «chaos-order». At α = 0.0321496, the alternation stops (Fig. 3,a). Up to α = 0.0321596, we see the modes of simple perio dic oscillations 14⋅2 0 .
In Fig. 5, we constructed the projections of a Poincaré map by the plane P = 0.2 for, respectively, the regular attractor of a simple periodic cycle 8⋅2 0 (α = 0.0321107) (a); the regular attractor of a quasiperiodic cycle ≈ 1⋅2 0 on a toroidal surface (α = 0.0321125) (b); the regular attractor математичне моДелювання біохімічних процесів The projections correspond to modes, whose type is determined by the calculated Lyapunov exponents. The location of the points of intersection of the plane by a trajectory corresponds to the topological configuration of the trajectory.
In Fig. 5,a, we show the intersection points for a 8-fold periodic cycle. These points are strongly repeated in the period of the given cycle. Figure 5,b presents the set of the points of intersection of the plane by a quasiperiodic cycle ≈ 1⋅2 0 . Their maximum density is observed in the central part, which is a region of attraction of a onefold cycle. In Fig. 5,c,d, we give, respectively, the Poincaré sections for modes in the case where the phase trajectory moves along the torus surface. In this case, the phase trajectory on the torus surface can be considered as a superposition of two motions: a rotation around «parallels» of the torus with frequency f 1 and a rotation around the «cylinder» forming the torus with frequency f 2 . The points of intersection of the plane P = 0.2 by the phase trajectory arise in regular time intervals equal to the period t = 1/f 1 . The shape of the Poincaré section depends on the ratio f 1 /f 2 . At α = 0.0321286, this ratio is irrational. Therefore, the phase trajectory is closed and covers the torus surface everywhere densely. The regular attractor of a quasiperiodic cycle ≈ n⋅2 0 is formed. The Poincaré section is a closed curve in the form of an ellipse (Fig. 5,с). At α = 0.032135, the ratio f 1 /f 2 = n 1 /n 2 is rational, and the Poincaré section consists of a finite number of points, which are located on a closed curve. In this case, the periodic motion with the period t = n 1 /f 1 = n 2 /f 2 occurs.
In other words, the phase trajectory is closed after n 1 rotations around «parallels» and n 2 rota-математичне моДелювання біохімічних процесів tions around «meridians» with the formation of a regular attractor 47⋅2 0 . The Poincaré section contains 47 points (Fig. 5,d).
According to Fig. 5,b,c, we constructed the Poincaré maps for the variable e 1 in Fig. 6,a,b. For a quasiperiodic cycle ≈ 1⋅2 0 , the Poincaré map is the set of points with their maximum localization in the central region. For the regular attractor of a quasiperiodic cycle ≈ n⋅2 0 , the Poincaré map is an oval characterizing the repetition of oscillations on the torus.
As an example for the visual representation of a complicated structure of the regular attractor 47⋅2 0 on a toroidal surface (α = 0.032135), we present its three-dimensional projections in the coordinates (e 1 ,ψ,G) (a) and (B,ψ,o 2 ) (b) in Fig. 7,a,b. Despite a high strokeness, the trajectory returns strictly regularly at any point of the cycle in the time interval strictly equal to the cycle period.
The attractors under consideration are presented in Table. The strange attractor 8⋅2 x a b c d (α = 0.03211295) is not unique in the given region. Strange attractors can be found at other values of α at a more thorough study. But all they have a common mechanism of formation of their structure. In accordance with the values of Lyapunov exponents (λ 3 -λ 10 ) < 0, the trajectories of the given attractor exponentially converge in the process of evolution.
In the same directions, the element of the phase volume shrinks. But since λ 1 > 0, the trajectories diverge in this specific direction. Moreover, the element of the phase volume is extended in this direction, though it shrinks during the evolution on the whole, Λ < 0. Such a motion results in that the phase trajectories form folds and mix with one another. In such situations, the smallest fluctuation will cause the randomness of a trajectory. The elongation and the folding of a chaotic attractor eliminate the initial information and replace it by a new one. Thus, the course of evolution of the system becomes unpredictable. The system can be characterized by the deterministic chaos as a model of chaos arising in the real system. The metabolic process in a cell adapts itself to these conditions, by preserving its stability in certain limits of the stability of a strange attractor. The further complete investigation of modes of the metabolic process running in cells within the given model and the consideration of conditions for the appearance of a deterministic chaos will be performed in the subsequent works.
With the use of the mathematical model of a cell arthrobacter globiformis, we have studied the internal dynamics of its metabolic process in the oscillatory mode. We have found the scenarios of formation and destruction of regular and strange attractors with various periods and types depending on the dissipation of the kinetic membrane potential. The limits of phase-parametric characteristics of the regions, where the bifurcations and the transitions «order-chaos», «chaos-order» and «orderorder» arise, are determined. By calculating the total spectra of Lyapunov indices, divergences, Poincaré sections, and Poincaré maps, we studied the regularity of the sequence of attractors appearing on a toroidal surface. In particular, we found a strange attractor, whose structure is formed with the help of a fold. In such a mode, the metabolic process in a cell adapts itself and conserves its stability characteristic of a strange attractor. математичне моДелювання біохімічних процесів