A mathematical model of the metabolism of a cell. Self-organization and chaos
UUDC 577.3 PACS 05.45.-a, 05.45.Pq, 05.65.+b
A mathematical model of the metabolism of a cell. Self-organization and chaos
V.I. Grytsay a,* , I.V. Musatenko b a Bogolyubov Institute for Theoretical Physics, Kyiv, Ukraine E-mail: [email protected] b Taras Shevchenko National University of Kyiv, Faculty of Cybernetics, Kyiv, Ukraine E-mail: [email protected]
Abstract:
Using the classical tools of nonlinear dynamics, we study the process of self-organization and the appearance of the chaos in the metabolic process in a cell with the help of a mathematical model of the transformation of steroids by a cell
Arthrobacter globiformis . We constructed the phase-parametric diagrams obtained under a variation of the dissipation of the kinetic membrane potential. The oscillatory modes obtained are classified as regular and strange attractors. We calculated the bifurcations, by which the self-organization and the chaos occur in the system, and the transitions “chaos-order”, “order-chaos”, “order-order,” and “chaos-chaos” arise. Feigenbaum’s scenarios and the intermittences are found. For some selected modes, the projections of the phase portraits of attractors, Poincaré sections, and Poincaré maps are constructed. The total spectra of Lyapunov indices for the modes under study are calculated. The structural stability of the attractors is demonstrated. A general scenario of the formation of regular and strange attractors in the given metabolic process in a cell is found. The physical nature of their appearance in the metabolic process is studied.
Keywords: mathematical model, metabolic process, self-organization, phase portrait, deterministic chaos, regular attractor, strange attractor, bifurcation, Poincaré section, Poincaré map, Lyapunov indices.
1. Introduction
In the present work, we continue the study of the mathematical model of the metabolic process in a cell
Arthrobacter globiformis . It is based on the process of transformation of steroids in a bioreactor, which is well investigated in experiments [1]. The constructed mathematical model allows us to determine the internal and external parameters, with which the model describes the stationary modes of a bioreactor. The studies within the model showed that autooscillations must appear in the biochemical reaction under certain conditions [2-17]. These autooscillations predicted as early as in 1985 [2] were found experimentally in [18, 19]. Analogous autooscillations are observed in the processes of photosynthesis, glycolysis, variations of the calcium concentration in a cell, oscillations in heart muscle, and other biochemical systems [20-24]. The study of such autooscillations will allow one to investigate the internal dynamics of metabolic processes in cells, to find the structural-functional connections in a cell, by which its vital activity runs, and to clarify the evolution of the formation of these connections. The application of the mathematical apparatus of nonlinear dynamics to the study of metabolic processes will allow one to develop the general methods of synergetics considering the physical laws of self-organization in the Nature.
2. Mathematical Model
The mathematical model of the metabolic process running in a cell
Arthrobacter globiformis at the transformation of steroids is constructed according to the general scheme of this process presented in Fig. 1. The model is based on the results of experimental studies of the process under flowing-hrough conditions with a fermenter in porious granules with immobilized cells
Arthrobacter globiformis [3, 4]. Fig. 1. General scheme of the metabolic process in a cell
Arthrobacter globiformis . The variation of the concentration of hydrocortisone ( G ) is described by the equation .3)()1(123 0 GGVEVlGN GdtdG ayg --++= (1) Under the action of the diffusion and the flow into pores of a macroporous granule to cells, hydrocortisone comes to the region of localization of the enzyme 3-ketosteroid- D -dehydrogenase ( E ) (term yg
23 0 ++ GN G ) and is transformed by this enzyme into prednisolone (term )()1(1
GVEVl ). A part of hydrocortisone is taken out from the biosystem by the flow (term G a ). Here and below, the function )( XV characterizes the adsorption of the enzyme in the region of local binding into active complexes; ).1/()( XXXV += The variation of the concentration of prednisolone ( P ): .4)()()2(2)()1(1 PPVNVEVlGVEVldtdP a--= (2) Prednisolone formed in the process (term )()1(1
GVEVl ) is transformed by the enzyme b
20 -oxysteroid-dehydrogenase ( E ) to its b
20 -oxyderivative (term )()()2(2
PVNVEVl ). Under the action of a flow (term P a ), a part of prednisolone goes out into the external solution. The variation of the concentration of b -oxyderivative of prednisolone ( B ): .5)()(1)()()2(2 BBVVkPVNVEVldtdB ay --= (3) The increase of the concentration of B occurs as a result of the transformation of prednisolone (term )()()2(2 PVNVEVl ). Its decrease is due to the use of b
20 -oxyderivative by cells in one of the possible modifications of the Krebs cycle (term )()(1
BVVk y ), which increases the level of HNAD (cid:215) . Under the action of a flow (term B a ), B is washed out into the external solution. The variation of the concentration of the oxidized form of 3-ketosteroid- D -dehydrogenase ( E ): -+++-+= )11(21 2101 mNPN mNPGGEdtdE b .11)()1(4)()1(1 EQVeVlGVEVl a-+- (4) The biosynthesis of the enzyme is described by the term )11(21 2 mNPN mNPGGE +++-+b , which is defined by the activation by the substrate G and the inhibition by the reaction products P and N . The decrease of the concentration of this form of the enzyme in the process of transformation of hydrocortisone is given by the term )()1(1 GVEVl , and its increase in the process of reduction of the respiratory chain corresponds to the term )()1(4
QVeVl . The inactivation of the enzyme due to the proteolysis is described by the term 11 E a . The variation of the concentration of the reduced form of 3-ketosteroid- D -dehydrogenase ( e ): .11)()1(1)()1(41 eGVEVlQVeVldtde a-+-= (5) Its level decreases in the process of reduction of the respiratory chain (term )()( QVeVl - ) and due to the inactivation (term e a ) and increases at the transformation of hydrocortisone (term )()1(1 GVEVl ). The variation of the level of the oxidized form of the respiratory chain ( Q ) ),()(71)()1(6)()1()2()2(6 NVQVlQVeVlVOVQlVdtdQ ---= y (6) where )21/(1)()1( yy += V . We accept that the concentration of menaquinone 200 =+ qQ , where q is the reduced form of the respiratory chain. The respiratory chain is oxidized by oxygen (term )()1()2()2(6 y VOVQlV - ) and is reduced with the help of 1 e (term )()1(6 QVeVl - ) and due to the high level of HNAD (cid:215) (term )()(7
NVQVl - ) . The variation of the concentration of oxygen ( O ): .27)()1()2()2(25 202 OVOVQlVON OdtdO ay ---+= (7) Under the action of a flow (terms
25 2 ON O + and 27 O a ), the level of aeration of a cell is changed. The concentration of oxygen decreases at the oxidation of the respiratory chain (term )()1()2()2( y VOVQlV -- ). The variation of the concentration of b -oxysteroid-dehydrogenase ( E ): -+-++= )21(22 2202 BN BNNPPEdtdE bb EPVNVEVl a-- (8) The increase of the level of the given enzyme occurs due to the biosynthesis: )21(22 220
BN BNNPPE +-++ bb . Prednisolone and
HNAD (cid:215) are activators of this process, and b
20 -oxyderivative is an inhibitor. The decrease of the level of the given enzyme occurs as a result of the inactivation ( E a- ) and the process of transformation of prednisolone ( )()()2(10 PVNVEVl - ). The variation of the concentration of HNAD (cid:215) ( N ) of a cell: +--= )()(7)()()2(2 NVQVlPVNVEVldtdN .64 010)(2
NNN NKBVk ayy -++++ (9) he level of the co-enzyme N decreases in the process of transformation BP ⇒ , in the process of reduction of the respiratory chain ( )()(7 NVQVl - ), and due to a flow ( N a - ). It increases at the use of B by cells in the Krebs cycle as a substrate ( yy + KBVk ) and in the presence of endogenous substrates (
NN N + ) in the environment. The variation of the level of kinetic membrane potential ( y ): ayy -+= )()(8)()1(5 QVNVlGVEVldtd . (10) The kinetic membrane potential arises at the transformation of hydrocortisone ( )()1(5
GVEVl ) and the reduction of the respiratory chain ( )()(8
QVNVl ) at a high level of HNAD (cid:215) and decreases due to other metabolic processes ( ay- ). The variation of the level of y changes its regulatory role (1), (3), (6), (7), (9). If the potential is high, the respiratory chain is blocked and held in the reduced state. The main parameters of the system, with which we fit the relevant experimental data, are as follows: ;2.011 === kll ;27.0102 == ll = l ; ;5.064 == ll ;2.17 = l ;4.28 = l ;5.12 = k ;310 = E ;21 = b ;03.01 = N ;5.2 = m = a ; ;007.01 = a = a ; ;2.120 = E ;01.0 = b ;12 = b ;03.0 = N ;02.02 = a ;019.00 = G ;23 = N ;2.02 = g ;014.05 = a ;001.07643 ==== aaaa ;015.020 = O ;1.05 = N = N ; ;14 = N = K . The study of solutions of the given mathematical model was carried out with the help of the theory of nonlinear differential equations [25-27]. In the numerical solution of this autonomous system of nonlinear differential equations, we used the Runge--Kutta--Merson method. The accuracy of calculations was set to be - . To attain the reliability of a solution, when the system passes from the initial transient phase onto the asymptotic solution with an attractor, the duration of calculations was taken to be . For this time interval, the trajectory “sticks” onto the appropriate attractor. The various types of autooscillatory modes are studied with the help of the construction of exact phase-parametric diagrams. We found the scenarios of appearance of bifurcations at the transition of the dynamical process from one type of an attractor to another one. For the most characteristic modes, we calculated the total spectra of Lyapunov indices (Table 1). To construct a phase-parametric diagram, we used the method of section. In the phase space of trajectories of the system, we place a cutting plane with P = 0.2. Such choice is explained by the symmetry of oscillations relative to this point of this variable in multiple modes. If the trajectory P ( t ) crosses this plane in a certain direction, we mark the value of chosen variable (e.g., G ) on the phase-parametric diagram. In such way, we have the point corresponding to the section of a trajectory by the two-dimensional plane. If the multiple periodic limiting cycle appears, we obtain a number of points, which will be coincide in a period. If a deterministic chaos arises, the points of intersection of trajectories by the plane will be placed chaotically. In order to uniquely identify the form of an attractor for the chosen points, we calculated the total spectrum of Lyapunov indices and determined their sum ∑ =L j j l (see Table 1). The calculation was carried out by Benettin’s algorithm with orthogonalization of the vectors of perturbation by the Gram--Schmidt method [26, 28, 29].
3. Results of Studies
We now consider the dynamics of modes within the mathematical model (1)-(10) under a variation of the dissipation of a kinetic membrane potential a (10) [16, 17]. We found the autooscillatory and chaotic modes with various multiplicities. The projections of their phase portraits have a characteristic form shown in Fig. 2,a,b. Fig. 2. Projections of the phase portraits of regular attractors; a – autoperiodic cycle (cid:215) for a = 0.033; b – quasiperiodic cycle (cid:215)» for a = ˛a (0.032159, 0.032166). Fig. 3. Bifurcation diagram of the system for ˛ a (0.032159, 0.32166). For ˛ a (0.0321590, 0.03215960), the regular attractor of the 14-fold period (cid:215) is kept in the system. For a = 0.03215961, we observe the appearance of the period doubling bifurcation with the generation of the regular attractor 1214 (cid:215) (Table 1). Then for a = 0.03215962, there arises the bifurcation of the generation of a two-dimensional torus (the Neimark bifurcation). The configuration of kinetic curves is instantly changed, and the quasiperiodic attractor with n -fold period is established on the toroidal surface (cid:215)» n ( t ) (Figs. 4,a and 5,a). a b Fig. 4. Kinetic curve of the variable )(1 te ; a - regular attractor of the quasiperiodic cycle » n*
2 on the toroidal surface for a =0.03215962. b - regular attractor (cid:215) for a =0.032162. Fig. 5. Projections of phase portraits; a – regular attractor of the quasiperiodic cycle (cid:215)» n on the toroidal surface for a =0.03215962; b – strange attractor x (cid:215) for a =0.032164. s a increases, the given attractor loses the stability, by passing periodically to the 14-fold limiting cycle ( a = 0.032160), which corresponds to the gaps in Fig. 3,a. In addition, other various multiple modes arise. For example, for a = 0.032161, 0.0321615, and 0.032162, the regular attractors (cid:215) , (cid:215) , and (cid:215) appear, respectively (Fig. 4,b). As a increases, we see the appearance of bifurcations of the limiting cycle. Moreover, the instant structural rearrangement of the type “order-order” occurs; i.e., as a result of the self-organization, the regular attractor of some form is replaced instantly by a regular attractor of some other form. In this case, the trajectories leave the region of attraction of the attractor and are drawn in the region of attraction of another regular attractor. The interesting scenario of the metabolic process is observed in the interval ˛ a (0.0321626, 0.032164). In Fig. 6, we present a magnified part of the bifurcation diagram in Fig. 3. Fig. 6. Phase-parametric diagram of the system for ˛ a (0.0321626, 0.32164), where Feigenbaum’s scenario is observed. At the beginning of the interval at a = 0.0321626, the regular attractor (cid:215) is formed on the toroidal surface. For j a = 0.03216276, the bifurcation yields the doubling of the period, and the regular attractor (cid:215) arises on the toroidal surface. For + j a = 0.03216346 and + j a = 0.03216361, we see the attractors (cid:215) and (cid:215) , respectively. This sequence of bifurcations satisfies the relation »+-+ -+¥fi jj jjt aa aa d = 4.669211660910… characterizing the infinite cascade of bifurcations at the transition to a deterministic chaos. Thus, as the coefficient of dissipation a increases in this region, the period of a complicated regular attractor on the torus is doubled by Feigenbaum’s scenario [37-40]. The further increase in a causes a deviation from the given scenario and the formation of the strange attractor x (cid:215) ( a =0.032164, Fig. 5,b) as a result of the intermittency. But then, for a =0.032174, the strange attractor x (cid:215) appears (Fig. 7,b). In the interval a = (0.032164, 0.032174) as a result of the intermittency of these chaotic cycles, we observe the transition between them: x (cid:215)« . In Fig. 7,a for a =0.032165, we show a projection of the phase portrait of a mutual transition of the given strange attractors. Figure 8 presents the kinetic curve for the variable (1 te for tis mode. We observe the transition “chaos-chaos”: x (cid:215)« . Moreover, the strange attractor x (cid:215) on the left and the strange attractor x (cid:215) on the right move toward each other. Since there are no other attractors of the system in this region, the trajectory is chaotically kept in the region of attraction of the strange attractor x (cid:215) or the strange attractor x (cid:215) Under the effect of bifurcations, the trajectory is aperiodically drawn in one of the regions of the given strange attractors after the transient process. According to the values of higher Lyapunov indices (Table 1), the formed limiting set is unstable by Lyapunov. Fig. 7. Projections of the phase portraits; a – strange attractor of the mutual transition x (cid:215)« for a =0.032165; b – strange attractor x for a =0.032174. Fig. 8. Kinetic curve of the variable )(1 te of the mutual transition of the strange attractors x (cid:215)« for a =0.032165. For the given strange attractor, we constructed a projection of the section by the plane P = 0.2 and the Poincaré map in Fig. 9,a,b. The choice of a cutting surface was made to attain the maximum number of intersections of the given component and the phase trajectory P(t) , as the former decreases, without contacts. Fig. 9. Projection of the section by the plane P = 0.2 (a) and Poincaré map (b) of the strange attractor formed during the mutual transition x (cid:215)« for a =0.032165. The obtained points of intersections and the Poincaré maps are grouped along several curves that form a geometric self-similarity. On the projection, we see clearly the fractality of this strange attractor. In addition, these curves do not create a quasistrip structure. Their number increases permanently with the duration of numerical integration of the system. This testifies to the impossibility of any reduction of the given complicated mathematical model to some one-dimensional discrete approximation without loss of the information about the dynamics of the metabolic process in a cell. We note that the general scheme (Fig. 2) includes only the main parts of the metabolic process running in any cell with substrate-enzyme reactions and in the respiratory chain. Therefore, the model gives a rather general qualitative representation of the dynamics of the internal self-organization of the metabolic process in a cell. Table 1. Total spectra of Lyapunov indices for attractors of the system under study ( l - 10 l are not important for our investigation). a Attractor l l l L (cid:215) .000056 -.000214 -.003250 -.898509 0.0321596 (cid:215) .000040 -.000142 -.003306 -.898550 0.03215961 (cid:215) .000078 -.000150 -.003394 -.899865 0.03215962 )(02 tn (cid:215)» .000063 .000026 -.000274 -.905553 0.032160 0214 (cid:215) .000040 -.000146 -.003365 -.899368 0.032161 (cid:215) .000051 -.000142 -.000123 -.905352 0.0321615 (cid:215) .000062 -.000596 -.000576 -.902277 0.032162 (cid:215) .000064 -.000171 -.000155 -.905320 0.0321626 )(027 t (cid:215) .000063 -.000097 -.001180 -.902078 0.03216276 )(127 t (cid:215) .000062 -.000005 -.001267 -.902189 0.03216346 )(227 t (cid:215) .000047 .000025 -.001252 -.902056 0.03216361 )(327 t (cid:215) .000048 -.000023 -.001265 -.902267 0.032164 x (cid:215) .000367 .000018 -.001641 -.902164 0.032165 x (cid:215)« .000363 -.000004 -.001598 -.904005 .032174 x (cid:215) .000693 .000020 -.003534 -.901422
4. Conclusions
We have constructed a mathematical model of the metabolic process in a cell
Arthrobacter globiformis at the transformation of steroids. With the help of the given model, we have found the autooscillations in agreement with experiment, which show the complicated internal dynamics in a cell. The model is optimized by the number of variables of the system required for a qualitative description of the metabolic process under study. The given model involves the general regularities characteristic of any cell consuming a substrate, on the whole. The autooscillations arise on the level of the substrate-enzyme interaction with participation of the redox process in the respiratory chain and characterize the times of such interactions. At the synchronization of the given processes, the autooscillations characterizing the self-organization of the metabolic process on the whole are observed. At the desynchronization of the given processes, we see the adaptation of the metabolic process in a cell to varying external conditions in the environment with conservation of its functionality. The scenario of the transitions “order-chaos”, “chaos-order”, “order-order”, and “chaos-chaos” is studied with the help of Poincaré sections and maps. The total spectra of Lyapunov indices are calculated, and the structural stability of the obtained attractors is studied. Feigenbaum’s scenario and the Neimark bifurcation are found. The results will allow one to carry on the search for metabolic oscillations in a cell and to clarify the physical laws of self-organization.
Acknowledgement
The work is supported by the project N 0112U000056 of the National Academy of Sciences of Ukraine.
References A. A. Akhrem and Yu. A. Titov.
Steroids and Microorganisms , Nauka, Moscow, 1970 (in Russian). 2.
V. P. Gachok and V. I. Grytsay . Kinetic model of macroporous granule with the regulation of biochemical processes.
Dokl. Akad. Nauk SSSR
V. P. Gachok, V. I. Grytsay, A. Yu. Arinbasarova, A. G. Medentsev, K. A. Koshcheyenko and V. K. Akimenko. Kinetic Model of Hydrocortisone 1-en Dehydrogenation by
Arthrobacter globiformis . Biotechn. Bioengin.
33: 661-667, 1989. 4.
V. P. Gachok, V. I. Grytsay, A. Yu. Arinbasarova, A. G. Medentsev, K. A. Koshcheyenko and V. K. Akimenko. A kinetic model for regulation of redox reactions in steroid transformation by
Arthrobacter globiformis cells.
Biotechn. Bioengin.
33: 668-680, 1989. 5.
V. I. Grytsay. Self-organization in the macroporous structure of the gel with immobilized cells. Kinetic model of bioselective membrane of biosensor.
Dopov. Nats. Akad. Nauk Ukr.
No 2: 175-179, 2000. 6.
V. I. Grytsay. Self-organization in a reaction-diffusion porous media.
Dopov. Nats. Akad. Nauk Ukr.
No. 3: 201-206, 2000. 7.
V. I. Grytsay. Ordered structure in a mathematical model biosensor.
Dopov. Nats. Akad. Nauk Ukr.
No. 11: 112–116, 2000. 8.
V. I. Grytsay. Self-organization of biochemical process of immobilized cells bioselective of membrane biosensor.
Ukr. J. Phys.
46: 124-127, 2001. 9.
V. V. Andreev and V. I. Grytsay. Modeling of inactive zones in porous granules of a catalyst and in a biosensor.
Matem. Modelir. , no. 2: 57-64, 2005. 10. V. V. Andreev and V. I. Grytsay. Influence of heterogeneity of diffusion-reaction process for the formation of structures in the porous medium.
Matem. Modelir.
17, no. 6: 3-12, 2005. 11.
V. I. Grytsay and V. V. Andreev. The role of diffusion in the active structures formation in porous reaction-diffusion media.
Matem. Modelir.
18, no. 12: 88-94, 2006. 12.
V. I. Grytsay. Unsteady Conditions in Porous Reaction-Diffusion.
Roman. J. Biophys.
17, no. 1: 55-62, 2007. 3.
V. I. Grytsay. The uncertainty in the evolution structure of reaction-diffusion medium bioreactor.
Biofiz. Visn.
No.2: 92-97, 2007 14.
V. I. Grytsay. Formation and stability of morphogenetic fields of immobilized cell in bioreactor.
Biofiz. Visn.
No. 2: 25-34, 2008. 15.
V. I. Grytsay. Structural instability of a biochemical process.
Ukr. J. Phys.
55: 599-606, 2010. 16.
V. I. Grytsay and I. V. Musatenko. Self-oscillatory dynamics of the metabolic process in a cell.
Ukr. Biochem. J.
85, no. 2: 93 – 104, 2013. 17.
V. I. Grytsay and I. V. Musatenko.
The structure of a chaos of strange attractors within a mathematical model of the metabolism of a cell.
Ukr. J. Phys.
58: 677-686, 2013. 18.
A. G. Dorofeev, M. V. Glagolev, T. F. Bondarenko and N. S. Panikov. Unusual growth kinetics of
Arthrobacter globiformis and its explanation.
Mikrobiol. , no. 1: 33-42, 1992. 19. A. S. Skichko and E. M. Koltsova. A mathematical model to describe the fluctuations biomass of bacteria.
Teor. Osnov. Khim. Tekhn.
40, no. 5: 540-550, 2006. 20.
E. E. Selkov. Self-oscillations in glycolysis.
Europ. J. Biochem.
4: 79-86, 1968. 21.
B. Hess and A. Boiteux. Oscillatory phenomena in biochemistry.
Ann. Rev. Biochem .
40: 237-258, 1971. 22.
A. Goldbeter and R. Lefer. Dissipative structures for an allosteric model. Application to glycolytic oscillations.
Biophys J.
12: 1302-1315, 1972. 23.
A. Godlbeter and R. Caplan. Oscillatory enzymes.
Ann. Rev. Biophys. Bioeng.
5: 449–476, 1976. 24.
Chaos in Chemical and Biochemical Systems . Field R., Györgyi L. (eds.), World Scientific, Singapore, 1993. 25.
V. S. Anishchenko.
Complex Oscillations in Simple Systems , Nauka, Moscow, 1990 (in Russian). 26.
S. P. Kuznetsov,
Dynamical Chaos , Nauka, Moscow, 2001 (in Russian). 27.
T. S. Krasnopol’skaya and A. Yu. Shvets,
Regular and Chaotic Dynamics of Systems with Bounded Excitation , R&C Dynamics, Moscow-Izhevsk, 2008 (in Russian). 28. I . Shimada and T. Nagashima.
Progr. Theor. Phys . 61: 1605-1616, 1979. 29.