Global Self-Organization of the Cellular Metabolic Structure

Background Over many years, it has been assumed that enzymes work either in an isolated way, or organized in small catalytic groups. Several studies performed using “metabolic networks models” are helping to understand the degree of functional complexity that characterizes enzymatic dynamic systems. In a previous work, we used “dissipative metabolic networks” (DMNs) to show that enzymes can present a self-organized global functional structure, in which several sets of enzymes are always in an active state, whereas the rest of molecular catalytic sets exhibit dynamics of on-off changing states. We suggested that this kind of global metabolic dynamics might be a genuine and universal functional configuration of the cellular metabolic structure, common to all living cells. Later, a different group has shown experimentally that this kind of functional structure does, indeed, exist in several microorganisms. Methodology/Principal Findings Here we have analyzed around 2.500.000 different DMNs in order to investigate the underlying mechanism of this dynamic global configuration. The numerical analyses that we have performed show that this global configuration is an emergent property inherent to the cellular metabolic dynamics. Concretely, we have found that the existence of a high number of enzymatic subsystems belonging to the DMNs is the fundamental element for the spontaneous emergence of a functional reactive structure characterized by a metabolic core formed by several sets of enzymes always in an active state. Conclusions/Significance This self-organized dynamic structure seems to be an intrinsic characteristic of metabolism, common to all living cellular organisms. To better understand cellular functionality, it will be crucial to structurally characterize these enzymatic self-organized global structures.


Introduction
One of the most important goals of contemporary biology is to understand the elemental principles governing metabolic structure as a whole which underlie the common design of all microorganisms and cells.
This global metabolic structure, conformed by the reactive interactions of thousands of biochemical species densely integrated through a labyrinthine web, represents one of the most complex dynamic systems in nature [1].
During the preceding two decades, different metabolic models have been studied intensively. Traditional models have focused on the kinetics of multi-enzyme systems by solving systems of differential equations and algebraic equations [2]. Petri's net theory, among others [3], has been applied to modelling metabolic pathways [4], decomposition of large metabolic networks into smaller subnetworks [5] and topological analysis of metabolic networks [6]. Large networks present many connections between the nodes, and their degree distributions follow a power law, so they can be considered as scale-free [7,8]. The presence of ''smallworld'' features [9] in scale-free networks is being studied [10,11]. Constraint-based modeling approaches, such as flux-balance analysis, has been applied in several metabolic networks [12,13]. Other mathematical models have been proposed to organize the networks both in its modular and hierarchical structure [14][15][16].
In an attempt to get a more accurate comprehension of the metabolic dynamic phenomena, we have proposed a dynamical system called ''dissipative metabolic networks'' or DMNs, which is basically formed by groups of enzymatic associations dissipatively structured and interconnected by fluxes and regulatory signals (allosteric and covalent).
Our model takes into account the fact that the cellular organization at the molecular level presents two relevant dynamic characteristics: the presence of enzymes aggregated in clusters and the emergence of dissipative catalytic patterns.
Experimental observations have shown that enzymes operate within metabolic pathways, may form functional catalytic associations and, thus, do not function in isolation of one another. Some of the first experimentally isolated enzymatic associations were, among others, the glycolytic subsystem [17], five enzymes from the cycle of the tricarboxylic acid [18], a triple multienzymatic-associate formed by the alpha-ketoglutarate dehydrogenase complex, the isocitrate dehydrogenase and the respiratory chain [19], and the complex formed by malate-dehydrogenase, fumarase and aspartate transferase [20]. Nowadays there are enough experimental data confirming the existence of numerous enzymatic associations belonging to metabolic routes, like lipid synthesis, glycolysis, protein synthesis, Krebs cycle, respiratory chain, purine synthesis, fatty acid oxidation, urea cycle, DNA and RNA synthesis, amino acid metabolism, cAMP degradation, etc. [21][22][23][24][25][26][27][28].
Association of various enzymes in large complexes (supramolecular organization) allows the direct transfer of their common intermediate metabolites (metabolic channelling). In addition, reversible interactions of enzymes with structural proteins and membranes are a common occurrence. This results in the existence of microcompartments within the soluble phases of cells. The microcompartmentation provides biophysical and biochemical mechanisms of physiological importance for the control of metabolic pathways [29][30][31][32].
The second consideration in our model is the presence of dissipative catalytic patterns. Each functional enzymatic association conforms a catalytic entity as a whole, in which spontaneously organized molecular oscillations may emerge.
In the far from equilibrium conditions prevailing inside the cell, the catalytic dynamics of enzymatic sets may present transitions between different stationary and oscillatory molecular patterns. When the enzymatic sets exhibit a rhythmic behavior, all the metabolic intermediaries oscillate with the same frequency but different amplitudes [33].
It is well known that these biochemical rhythms constitute a new type of supramolecular organization that may emerge in open systems far from equilibrium and was called dissipative structure by Prigogine [34].
It has been found experimentally that many enzymatic associations may cause oscillatory processes. The first systematic classification of cellular oscillatory behaviors gathered about 450 different kinds of biochemical rhythms [35], most of them corresponding to periodic oscillations [33,36].
Due to the importance of cellular enzymatic rhythms, we have taken into account this aspect in our model of DMNs, so that the simulated catalytic processes can present both stationary and oscillatory activity regimes.
In the dissipatively structured enzymatic associations, the existence of allosteric enzymes permits the interconnection among them. The catalytic activity of the allosteric enzymes is modulated through the noncovalent binding of a specific metabolite at a place of the protein different from the catalytic site, provoking alterations of the metabolic state in an interval of seconds. Such types of modulation may be both positive (activation of their catalytic rates) and negative (inhibitory modulators). The regulation by means of the covalent interactions can originate ''all-ornothing'' type answers [62].
In agreement with all these considerations, our DMNs model is composed of a set of catalytic elements (each of them represents a dissipatively structured enzymatic association called metabolic subsystem), which are connected by substrate fluxes and regulatory signals (allosteric and covalent modulations). The metabolic subsystems may present oscillatory and stationary activity patterns.
In a previous work [63], we have shown the rich variety of selforganised temporal patterns and global configurations that appears in the DMNs. One of them is the emergency of a functional configuration (similar to in vivo) in which a set of metabolic subsystems are always locked into active states (in eukaryotic cells they are mainly the Krebs cycle subsystem, the pyruvate dehydrogenase complex and the oxidative phosphorylation), whereas the rest of metabolic subsystems present dynamics of on-off changing states (glycolysis, beta-oxidation of fatty acids, amino acid degradations, gluconeogenesis, etc).
This type of functional metabolic structure has recently been shown by Almaas and colleagues for Escherichia coli, Helicobacter pylori, and Saccharomyces cerevisiae [12,64].
In this paper, we have analyzed around 2.500.000 different metabolic networks in order to investigate the conditions in which this global configuration emerges. We have found that this functional metabolic structure is an emergent property of all dissipative metabolic networks with a high number of enzymatic subsystems.

Model
Experimental observations have shown that the enzymes may form functional catalytic associations in which a new type of supramolecular self-organization may emerge (for more details, see the introduction section). These catalytic associations that operate within far from equilibrium conditions were called dissipative structures by Prigogine [34] and they conform a catalytic entity as a whole, in which the catalytic activity is autonomous with respect to the other enzymatic associations and molecular oscillations and steady states may emerge spontaneously. We have called metabolic subsystem to these groups of enzymatic associations dissipatively structured.
''Dissipative metabolic networks'' (DMNs) are dynamical systems basically formed by a given number of groups of dissipatively structured enzymatic associations (called metabolic subsystems or MSbs) interconnected by fluxes of substrates and regulatory signals (both allosteric and covalent).
The metabolic subsystems receive distinct input fluxes of substrates and three possible types of regulatory signals: activatory, inhibitory (graduated) and total inhibitory (all-or nothing type). Each MSb transforms the input fluxes and regulatory signals into one or several output activities.
In agreement with experimental observations, the output activity of the MSbs may be oscillatory or steady state [33] and comprise an infinite number of distinct activity regimes.
In our model, when the set of dissipatively structured enzymes shows an activity with rhythmic behavior the output activities present nonlinear oscillations with different levels of complexity, as could be expected in the cellular conditions ''in vivo''.
In order to be able to effect these different oscillatory or steady state patterns in each metabolic subsystem, the input-output conversion is made in two stages. First, the input fluxes are transformed into an intermediary activity by means of flux integration functions. In a second phase, the ''intermediary activity'' is modified by means of the ''regulatory signals integration'', which depends on the combination of the received regulatory signals. Each regulatory signal has an associated regulatory coefficient which defines the intensity of its influence.
The parameters q i,j , called ''influence coefficients'', that will be described in the materials and methods section, represent the influence of activatory and inhibitory modulators on the metabolic subsystems. These interactions correspond with the regulatory activity of the allosteric enzymes.
If the signal is of total inhibition (all-or nothing type), then we will use a parameter d called threshold value (the level of the enzymatic covalent regulatory activity). When a given threshold value is reached it inhibits completely the activity of the MSb.
On the other hand, biological processes with long-term correlations are of notable interest in the study of complex dynamics.
In a previous work with DMNs [63] we have analyzed different transitions generated by several metabolic networks and we have observed that these transitions exhibit long-term correlations (in the studied metabolic subsystems, their activity values depend to some extent on the previous ones).
In particular, we have characterized the behavior of several large complex transitions series generated in metabolic subsystems by means of the Hurst exponent H, which has been shown to be a robust and reliable test to detect the presence of long-term correlations [63,74].
For a random series with independent increments, it can be shown that H should be equal to 0.5. If H?0.5 this is indicative of the presence of ''persistence'' in the data, which means that the present state of the system is affected by previous states.
The studied time series generated by DMNs demonstrate complex transitions between the activities of the metabolic subsystems and we have found values for the Hurst exponent of 0.07,H,0.4, indicative of long term persistence over all the studied ranges. Values of H,0.5 are interpreted as characteristic of ''trend-reversing'' or ''antipersistence'' (a type of ''persistence''). The behavior of the time series tends to reverse itself, for instance, a decreasing trend in the past usually implies an increasing trend in the future and conversely, an increase in the past is likely to be followed by a decrease. The high reliability of these results was tested with exhaustive Monte Carlo simulations. All the series studied presented persistence and these results showed clearly that the complex transitions in the DMNs exhibit long-term memory phenomena.
In an attempt to study a possible influence of the past activity on the enzymatic self-organized global structures we have considered the b parameter (the past influence coefficient) in some of the DMNs, in the sense that the present activity of each metabolic subsystem may be affected by its previous activities (see more details in the materials and methods section).
When the enzymatic activity is considered at the molecular level, the model must take into account the parameters for this level of organization. Thus, for example, we have studied the selforganization of certain enzymes at the molecular level; concretely, we have modelled the glycolytic subsystem by means of a system of differential equations with delay. In these studies, we have taken into account different molecular control parameters such as the Michaelis constant, dissociations constants, non-exclusive binding coefficients, equilibrium constants between conformational states, etc. [75][76][77].
On the other hand, experimentally, the dynamic structure of the cellular metabolism as a whole seems to be characterized by presenting a functional global configuration in which a metabolic core formed by several sets of enzymes are always in an active state, whereas the rest of the molecular catalytic sets exhibit dynamics of on-off changing states [12,64].
Our main goal is to understand the cause for the emergence of this global supramolecular dynamic organization of enzymes. Since this main goal focuses on a superior level of organization than the molecular level, we will consider in this paper the emergent dynamic behaviors in the DMNs obtained by means of changes in: the level of enzymatic covalent regulatory activity, the level of allosteric activity, the number of metabolic subsystems, the flux topology, the topology of the regulatory signals and the flux values.

Dissipative metabolic networks dynamics
Among the different dynamic behaviors that the DMNs may exhibit we have considered three aspects for our analyses.
1-Local activity developed by each metabolic subsystem (periodic and stationary patterns). 2-Dynamic transitions between steady states and periodic behaviors. 3-Global configurations that the DMNs may adopt: All the subsystems are always in an on state. II. All the metabolic subsystems always present cycles of activity-inactivity. III. A certain number of subsystems are always locked in an on state conforming a metabolic core, while the rest of the subsystems are in an on-off switching state.

Study of networks with two metabolic subsystems
In our first study, we have considered the simplest situation, which corresponds to a metabolic network formed only by two subsystems (concretely, we have analyzed the network described in the example of the materials and methods section and figure 1). In this study (see Table 1) the d threshold value in the regulatory signals of total inhibition (the level of the enzymatic covalent regulatory activity) has been fixed as control parameter and we have taken the past influence coefficient b to be 0.
At small threshold values, for 0#d#0.12 the MSb1 presents a single oscillatory behavior of one-period. In the interval 0.12,d#0.836 the first subsystem is always active and presents different cycles of transitions with 2, 3, 6, 12… and more of 100 Finally, for 0.836,d#1 deterministic chaotic transitions can be observed. In this range of d values, the network formed by only two metabolic subsystems spontaneously auto-organizes, provoking the emergence of a very complex chaotic behavior in which each subsystem presents infinite transitions between different periodic patterns. In this situation, the MSb1 modifies uninterruptedly its activity so that it never repeats itself for arbitrary long periods of time.
The MSb2 for 0#d#0.17 is inactive (small threshold values d represent high covalent regulatory activity). For 0.17,d#0.67 the second subsystem presents cycles of activity-inactivity with different patterns of transitions between steady states and periodic behaviors. For 0.68#d#1 the MSb2 is locked in an active state. The same as with the first subsystem, for 0.837#d#1 deterministic chaotic transitions between steady states and periodic behaviors can be observed, that is, the MSb2 spontaneously exhibits infinite transitions between different oscillatory periodic behaviors and steady states.
In figure 2, some dynamical transitions corresponding to the mean amplitude A 0 and the corresponding periodic and stationary patterns developed by the two subsystems are shown.
The mechanism that determines these behaviors in both subsystems is not prefixed in any of the parts of the system. There is no feedback with oscillatory properties nor any other rules that determine the system to present complex transitions in the output activities of the metabolic subsystems. The dynamic behaviors which emerge spontaneously in the network have their origin in the regulatory structure of the feedback loops, and in the nonlinearity of the constitutive equations of the system (see materials and methods section for more details).
The introduction in our analyses of a new parameter, the b past influence coefficient, allows us to observe some changes in the dynamic behaviors of the metabolic nets formed by only two subsystems (see Table 2).
At small past influence coefficient b (b = 0.1), the threshold values d do not provoke qualitative changes in the network. So, for 0#d#1 the first subsystem is always on and chaos emerges between 0.8,d#1. In these conditions, the second subsystem exhibits three different activity states: for 0#d#0.2 the MSb2 is inactive, for 0.2,d#0.6 its state is on-off and for 0.6,d#1 the second subsystem presents an on state. Chaos emerges when 0.8,d#1.
When b = 0.2 the MSb1 is always on and chaos emerges in two parametric regions: for d = 0.3 and for 0.7,d#1. The second subsystem also presents, as the MSb1, chaotic behaviors under these conditions, however, for d = 0.3 the MSb2 exhibits an on-off state, and when 0.7#d#1 its state is on.  As b increases, the regions of chaotic behavior augment in correlation with d. For instance, when b = 0.5 and either 0.1,d#0.6 or 0.8,d#1 chaos emerge in the metabolic net, the first subsystem is always on and MSb2 is on-off in the first range and is on in the second one.
If b = 1 the net presents only one chaotic parametric region. The first subsystem is always on, and chaotic transitions emerges when 0.2,d#0.8. The second subsystem exhibits two different behaviors: for 0#d#0.2 and 0.7#d#1 its state is off, and when 0.2,d#0.7 its state is on-off; as occurs in the MSb1, chaotic transitions between steady state and periodic behaviors emerges for 0.2,d#0.7.
If we take d as the first reference parameter, it can be observed in Table 2

Study of metabolic networks with twelve metabolic subsystems
More interesting situations appear in the DMNs with more than one flux and one regulatory signal associated to each subsystem. In figure 3 we show a type of DMN formed by twelve subsystems (in the graphics only the interconnections by fluxes are reflected; the corresponding integration function parameters and the coefficient values of the regulatory signals are shown in Table 3).
Each subsystem may present three states: always on, always off and an on-off changing dynamic. If we take into account the dynamics followed by these subsystems, interesting functional configurations in the whole net can be observed.
Thus, when 0#d#0.29 all the subsystems are inactive. For 0.29,d#0.39, all the subsystems are locked into an on-off changeable state (figure 3a). In the network a qualitative change in the dynamical structure emerges for 0.39,d,0.68: various dissipative catalytic elements develop a set of subsystems locked in an on state, while the rest of metabolic subsystems are in an on-off changeable state (figure 3b). The modification of the control parameter in the range of 0.68#d#1, leads to a new phase transition in which all the subsystems are in an on state (figure 3c).
In Table 4 are represented different global configurations in DMNs formed by twelve subsystems in function of b and d. It is shown in bold the set of nets that present a self-organized global functional structure in which several subsets of enzymes are always in an active state (metabolic core) whereas the rest of molecular catalytic sets exhibit dynamics of on-off changing states. For very high values of b (0.9#b#1) never emerge DMNs characterized by presenting a metabolic core.    It can be observed also that in this kind of networks, for high covalent regulatory activity (0.1#d#0.4) nets with all metabolic subsystems in an on-off changing dynamic emerge (small threshold values d represent high covalent regulatory activity). Last, the DMNs with low levels of b and high values of d present all their metabolic subsystems always active (on).
In order to better understand the influence of the control parameters d and b on the global self-organizations in the DMNs, we have made a first study of different nets constituted by twelve metabolic subsystems with three regulatory signals and two input fluxes by subsystem. In this study all networks present the same flux configuration described in figure 3. The rest of parameters that configure the metabolic net, such as flux integration parameter, influence coefficients, and the regulatory interconnections were randomly chosen assigning the same probability to each type of regulatory signal.
A total of 1.210.000 different metabolic nets were built, the results of this study are shown in figure 4. In the calculations, 10,000 different metabolic networks were taken into account for each represented point (corresponding to determinate values of d and b), with 300 iterations per net. In these networks, for a given set of subsystems, random regulatory signals were settled. Each one of the regulatory feedbacks that connect the subsystems were randomly decided, and it was also decided randomly the type of regulation (activation, inhibition or total inhibition) as well as the parameters associated to each regulatory signal (see the materials and methods section for more details). The set of subsystems from whom a given subsystem receives fluxes was also decided randomly, (being thus established randomly the flux architecture of the net) as well as the parameters associated to the flux integration functions. The probability distribution that we used in the random selections was the uniform distribution. The control parameters were varied in steps of 0.1 in the unit interval [0, 1] and the criterion followed to determinate if a metabolic subsystem was on, off or changing on-off was that the subsystem had the corresponding state between the iterations 200 and 300.
In figure 4a, the percentage of DMNs with all their subsystems unable to change the state is shown (each subsystem is always on or is always off and never is in an on-off changeable state). It can be observed that if b is maximum (b = 1) all the percentages are greater than 50% (between 51.11% for d = 0 and 70.66% for d = 0.9). When the covalent regulation is very low (0.9,d#1) and 0,b#0.6, the biggest percentages are reached (between 64% and 97.6%). For 0.3,d#0.7 and 0.6,b#0.9 the lowest percentages are reached (all them are less than 10%).
The maximum percentage of dissipative networks with all its metabolic subsystems unable to change the state (97.6%) is obtained when d = 1 and b = 0.1, and the minimum (1.39%) is attained for d = 0.6 and b = 0.7.
In figure 4b the percentage of networks with all their subsystems in an off state is represented. They constitute a particular case of the DMNs showed in figure 4a and they are functionally unviable metabolic nets. When the past influence coefficient is maximum (b = 1) the biggest percentages are reached. These values range from 17.16% for d = 0.1 to 19.97% for d = 0.9. If 0.7#d#1 and b#0.7 all the percentages are less than 1%.
In figure 4c, the percentage of metabolic networks that present all their elements in an on-off regime is shown. The biggest percentages are obtained when b = 0.9 and 0.3#d#0.9. The lowest percentages (always less than 1%) are found for a value of d Table 4. Global functional configurations in the nets formed by twelve subsystems.  Finally, figure 4d shows the percentage of nets characterized by having a subset of dissipative metabolic elements locked into an active state (metabolic core) while the rest exhibit an on-off changeable state. If the covalent regulatory activity is minimum (d = 1), the percentage of nets presenting these functional configuration varies strongly as a function of b, in fact, the absolute maximum is 87.33% (b = 0.9) and the absolute minimum is 2.39% (b = 0.1). For a maximum covalent regulation (d = 0) the biggest percentage is 71.37% (b = 0.8) and the minimum percentage is 37.09% (b = 0). When 0,d,1, the highest percentages (about 80%) are attained for a given value of b between b = 0.5 and b = 0.7 and the minimum values correspond to b = 1. The 64.4% of the pairs (b, d) present a percentage over 50% for these kind of metabolic nets.
A total of 33.26% of the studied DMNs correspond to the case in which all its subsystems are unable to change the state, the 14.09% of the nets present all their elements in an on-off regime, and the 52.56% are metabolic nets characterized by having a subset of dissipative metabolic elements locked into an active state while the rest exhibit an on-off changeable state. The percentage of networks with all their subsystems in an off state (functionally unviable metabolic nets) is 4.95% (this percentage constitutes part of the 33.26% of the nets unable to change its state).
It can be observed that the global configuration characterized by having a metabolic core is the one that presents the maximum percentage when compared with the other two active global configurations.
In metabolic networks formed by twelve subsystems, complex dynamical transitions in the activities of each MSb are common and we have observed that chaotic transitions between steady states and periodic behaviors may emerge both in subsystems always on (for example, when d = 0.7 and b = 0) and in metabolic subsystems that exhibit an on-off changeable state (for example, when d = 0.8 and b = 0.3). Percentage of DMN that present, respectively, (a) all their subsystems unable to change the state (each subsystem is always on or is always off and never is in an on-off changeable state), (b) all subsystems in an off state (they constitute a particular case of the dissipative networks showed in figure 6a and they are functionally unviable metabolic nets), (c) all their elements in an on-off regime and (d) a subset of dissipative metabolic subsystems locked into an active state while the rest exhibit an on-off changeable state. In the horizontal axes the d threshold value (the level of the covalent regulatory activity) and the b past influence coefficient are displayed. In total, 1.210.000 different randomly constructed metabolic nets with 12 subsystems and two input flux by subsystem were studied. doi:10.1371/journal.pone.0003100.g004

Chaos in dissipative metabolic networks
Under experimental conditions, most of the local dynamic behaviors of the metabolic subsystems present periodical oscillatory patterns or steady states [33]. However, under cellular conditions it has been observed that some subsystems can present chaotic local activities.
In figure 5, two experimental chaotic behaviors are shown. The first one (figure 5a) corresponds to oscillations in citric acid cycle, a metabolic subsystem always active under cellular conditions [78] and the second (figure 5b) represents experimental calcium oscillations in Xenopus laevis oocyte carried out by our group. These chaotic oscillations correspond to a subsystem in an on-off changeable regime.
The blood serum from many animals contains a factor which activates a membrane receptor that is coupled to the phosphatidylinositol second messenger system and produces oscillatory currents. These currents are elicited by activation of Cl 2 channels sensitive to the intracellular Ca z 2 concentration [79]. The enzymatic activity bound to the membrane of the oocyte after the external stimulus of the Fetal Bovine Serum causes the dynamic mobilization of the intracellular calcium. This activity ceases soon after the exposure of the cell to the external agent.

Study of metabolic networks formed by high numbers of metabolic subsystems
Finally, in order to study the behavior of the DMNs depending on d (covalent regulation level) and n (number of subsystems), similar statistics were also carried out. In this case, networks formed by subsystems with one regulatory signal for each metabolic subsystem and one input flux were considered; likewise the same probability for each regulatory signal was assigned. Each net subsystem can have an outer input flux with a probability of 0.1, so that a metabolic network could have more than one outer input flux but at most one for each subsystem. Again, all parameters, the topology of the regulatory signals and flux interconnections (flux topology) were randomly configured. The b parameter was always taken equal to 0.3. For this study (figure 7) 1.250.000 different metabolic networks were analyzed.
In figure 7a, the percentage of metabolic nets with all their subsystems unable to change their state is shown. In this kind of nets, each subsystem is always on or is always off, and never is in an on-off switching state. It can be observed that for a fixed d there is a remarkable trend towards a decrease in this percentage when the number of subsystems that conforms the net is increased. When n = 2, the mean percentage with respect to d is 75.5% and for n = 50, it is 29.6%. With respect to the behavior of percentages when the number n of subsystems is fixed, two different clusters are observed. For n up to 35 subsystems, the percentage decreases slightly with d when d is less than 0.8, and increases starting from this value. In this cluster, the maximum percentage (83.31%) is obtained when d = 0.1 and n = 2, and the minimum percentage (24.81%) is obtained when d = 0.7 and n = 35.
For a number of metabolic subsystems greater than 35, the percentage decreases steadily with respect to d; this decreasing is very remarked for d.0.5 (that is, when the covalent regulation is low, the percentage of nets with all their subsystems unable to change their state diminishes notably). In this second cluster, the maximum percentage (41.92%) is obtained when d = 0.1 and n = 40, and the minimum percentage (21.65%) is obtained when d = 0.9 and n = 50.
The DMNs in which all subsystems are off constitute a particular case of the nets in which all the metabolic subsystems are unable to change their state (figure 7b). It is observed that, when d is fixed, the percentage goes asymptotically to zero. For a given value of n, the percentage decreases slowly with the increment of the d value. The maximum percentage (56.37%) is obtained when n = 2, and its minimum (0.32%) is reached when n = 50.
In figure 7c the percentage of metabolic nets with all its dissipative elements in an on-off switching state is shown. For d fixed, the percentage increases when the number of subsystems is between 2 and 4, and then it decreases, tending to 0 asymptotically with the number of subsystems. The global maximum (30.9%) is reached for d = 0.9 and n = 3, and the global minimum (0.11%) is reached for d = 0.1 and n = 50.
In figure 7d the percentage of DMNs that develop a metabolic core of subsystems locked into an active state while the rest of metabolic subsystems exhibit an on-off changeable state is shown. When d is fixed, the percentage grows pronouncedly with the number of subsystems. The rank of growth (that is, the difference of the percentage between n = 50 and n = 2) is 59.3 for d = 0.1 and 76.3 for d = 0.9. For a given n with n.35 the percentage increases notably with d.
The minimum percentage (0.61%) is attained when d = 0.1 and n = 2, and the maximum (77.71%) is obtained when d = 0.9 and n = 50.
It can be observed that an increase in the number of subsystems that constitute the net results in a decrease in the percentage of nets functionally non-viable (with all their subsystems off) as well as the nets with all their subsystems unable to change state and the nets with all their subsystems changing. Nevertheless, this increment in the number of subsystems suggests an asymptotic trend to reach a percentage of 100% of the nets presenting a global configuration characterized by having a metabolic core of subsystems locked into an active state while the rest of metabolic subsystems exhibit an on-off switching state.
These data seem to indicate that the fundamental element for the emergence of a global functional configuration characterized by presenting a metabolic core is a high number of metabolic subsystems.

Discussion
In order to study the emergent behaviors in metabolic structures formed by dissipative enzymatic associations (metabolic subsystems) connected by substrate fluxes and regulatory signals (allosteric and covalent interactions) we have used a type of dynamic system, which we call ''dissipative metabolic networks'' or DMNs.
Three types of basic emergent behaviors can be distinguished:

I.
Dynamic transitions corresponding to the mean amplitude, amplitude and frequency in the activities of each subsystem. II. Activity patterns developed by the metabolic subsystems. III. Global functional configurations in the network.
Upon the numerical analysis, a number of different qualitative types of transitions among the activity patterns of the subsystems can be observed: steady state-steady state, steady state-periodic regime, changes between different periodic regimes and chaotic transitions.
Periodic and chaotic transitions in these activity patterns of the subsystems are common. For the case of periodic cycles of transitions, the subsystems run repeatedly through the same set of states, resulting in a cycle of distinct periodic oscillations and steady states. Chaotic behaviors can be observed even in very simple nets of two or more subsystems with a single flux and a unique regulatory signal sent by a catalytic element in the network.
Each metabolic subsystem may present one of the three general states: always on, always off, or an on-off changing dynamics. If we consider the dynamic behaviors of a metabolic network regarding only these three general states, interesting global functional configurations in the overall of the net can be observed.
In fact, for a specific parameter value (for example, the threshold value in the regulatory signals of total inhibition, i.e. the level of enzymatic covalent activity) all the subsystems can be in an on state. In this situation, the dynamics are restricted to the possible changes in the variables of the active subsystems.
Further modifications in the control parameter may lead to a new global configuration in which all their dissipative subsystems are in an on-off changeable state.
A new variation in any control parameter may provoke a qualitative change in the metabolic net, resulting in the emergence of a new global configuration, in which a certain number of subsystems are locked in an active state (metabolic core) while the rest of the subsystems remain in an on-off changing dynamics. Besides, all these metabolic subsystems, both those locked in an always active state and those which are permanently in a process of activation and inhibition, exhibit local dynamics with transitions between different steady states and oscillatory behaviors.
In the prevailing conditions inside the cell, catalytic dynamics seem to show a similar structure to this kind of metabolic global configuration. The cellular metabolism presents a metabolic core formed by a set of catalytic associations always in active states (in eukaryotic cells they are mainly the tricarboxylic acid cycle, the pyruvate dehydrogenase complex and the oxidative phosphorylation) while the rest of the catalytic subsystems are in an on-off changing state (b-oxidation of fatty acids, amino acid degradations, glycolysis, gluconeogenesis, etc).
Since this global metabolic dynamic may be a genuine and universal functional configuration of the cellular metabolic structure common to all living cells, we have focused our efforts on analyzing the elements which may determine its emergence.
With this purpose, in a first study, we have analyzed 1.210.000 different randomly constructed metabolic nets with only 12 subsystems and with three regulatory signals and two input fluxes subsystem.
The numerical results show that more than 50% of the analyzed nets present a global configuration characterized by exhibiting a metabolic core.
Regarding the rest of the global metabolic configurations, apart from appearing in smaller percentages than the previous one, it is observed that these percentages diminish even more when the level of covalent regulation is not very high, as could be expected in living cell conditions (here small threshold values d represent high covalent regulatory activity).
In fact, when the nets with all their subsystems unable to change its state are studied (this happens for the 33.26% of the networks), it is noted that the smallest number of these nets is obtained when the level of covalent regulations is intermediate and the influence of the past is high (but not maximum, that is b,1).
In the particular case of nets with all their subsystems off (that case holds for the 4.95% of the nets, and constitutes part of the 33.26% of the nets unable to change its state), it can also be observed that the percentage of such nets is less than 1% when the level of covalent regulation is intermediate or low and the influence of the past is not very high (b,0.7).
As in the previous cases, when nets in which all the subsystems are in an on-off changeable state are considered (this happens for the 14.09% of the nets), it is observed that the lowest percentages are obtained when the level of covalent regulation is not high and the influence of the past is low.
The nets characterized by having a metabolic core, appear in the highest percentages, which become particularly high (up to 80%) when the level of covalent regulations is not high and b is low.
In the light of these results, despite the fact that the control parameters b and d may cause important changes in the dynamic behaviors of the nets, it seems clear that none of these control parameters are determinant in the emergence of a global metabolic configuration characterized by a set of subsystems locked into an active regime while the rest is in an on-off switching state in any generic metabolic net.
In order to better understand which elements could determine the emergence of this global configuration, we have carried out simulations including 1.250.000 randomly constructed nets, using the number n of subsystems that conform each network and the threshold value d (the enzymatic covalent regulatory activity level) as control parameters.
The obtained results seem to show that the percentage of nets with all their subsystems unable to change the state decreases remarkably with increased numbers of subsystems. This decrement is more notable when the level of covalent regulation in these nets is not high.
In the particular case of nets with all their subsystems in an off state (functionally unviable metabolic nets) their percentage tends asymptotically to zero. In fact, when the number of subsystems is 50, the percentage is negligible.
An analysis of the nets with all their subsystems in an on-off changeable state shows clearly that this percentage tends also to zero.
Finally, the study of the remaining nets (the ones with a metabolic core) shows a rapid increase in their percentage as a function of the number of subsystems. These values are particularly high when the level of covalent regulation is moderate. These results suggest that there is an asymptotic trend leading towards 100% in the percentages of this global functional configuration when the number of subsystems increases.
In conclusion, the fundamental element for the emergence of a global metabolic configuration characterized by presenting a metabolic core is the number of subsystems. It seems that this global configuration is an emergent property of all the DMNs with a high number of subsystems.
It is possible that the complexity inherent to a biochemical system conformed by a high number of dissipative subsystems provokes the spontaneous emergence of this fundamental aspect of the cellular metabolic structure.
Unicellular organisms present a kind of metabolic cellular structure characterized by having a dynamic functional configuration in which a small group of metabolic subsystems are always active whereas the rest present on-off catalytic dynamics. Our numeric analysis seems to show that this global dynamic selforganization of the metabolism emerges spontaneously in the cellular metabolic space, as a consequence of the fact that the cells possess catalytic structures conformed by a big number of metabolic subsystems.
All living cellular organisms seem to display some metabolic subsystems that are always active, (these are involved in the permanent production of ATP), whereas the rest of the subsystems present on-off reactive dynamics. Taking all this into account, we can conclude that this global dynamic organization of the reactive cellular processes is possibly a genuine characteristic of the metabolism, common to all living cellular organisms, which could be basic and fundamental in the regulation of the most elementary cellular processes.
By means of modelling approaches, such as the constraint-based analysis applied to metabolic networks, E. Almaas, A.L. Barabási and their group of researchers, [12,64] have shown that the global organization of metabolic fluxes in certain bacterial cells is characterized by displaying a metabolic core of catalytic reactions always active under different growth conditions, embedded in a connected set of enzymatic reactions where some enzymatic pathways are eventually turned off completely.
These studies have been made in Escherichia coli, Helicobacter pylori, and Saccharomyces cerevisiae metabolism, showing also that most current antibiotics may interfere with the metabolic core [64]. The authors suggest that this global organization of the cellular metabolism ''probably represents a universal feature of metabolic activity in all cells, with potential implications for metabolic engineering.' ' The global configuration that emerges in the DMNs, characterized by displaying a metabolic core formed by a set of subsystems always active embedded in a set of subsystems in an on-off changeable state, seems to agree remarkably with the observations of E. Almaas, A.L. Barabási and their group of researchers. Our numeric results seem to show that this metabolic global organization is originated by the dissipative dynamics that spontaneously emerge in the cellular reactive structure formed by a network with a high number of enzymatic processes.
On the other hand, the emergence of frozen cores in dynamic networks has also been pointed out in Random Boolean networks. This interesting kind of dynamical systems was introduced in 1969 by S. Kauffman [81][82][83].
Under experimental conditions, most of the local dynamic behaviors of the metabolic subsystems exhibit periodical oscillatory patterns or steady states. However, it has also been observed that some subsystems can present chaotic activities. For example, that occurs in oscillations in citric acid cycle (an always active metabolic subsystem) and in calcium oscillations in oocytes where the membrane receptors are activated.
Our results also show that the mean amplitude, amplitude and frequency in the activities of certain subsystems may present chaotic patterns. These dynamic behaviors may emerge in the simplest nets with only two subsystems and in the most complex nets; likewise, in subsystems always on as well as in metabolic subsystems that exhibit an on-off switching state chaos can be observed.
The combination of chaotic dynamics in some subsystems with a stable functional configuration which presents a metabolic core may be of biological interest. The stability of this global configuration is necessary to ensure the maintenance of its metabolic structure. The existence of chaotic dynamics in the transitions of the subsystems activity may constitute an advantage in the self-regulatory control of the system, due to the sensitivity of chaotic behaviors to initial conditions; these dissipative metabolic networks with local chaotic patterns may permit fast responses during the cellular adaptation and in the regulation against perturbances.
This conception of the cellular metabolic structure as a complex dissipative catalytic network endowed with a stable global configuration in which a core of metabolic subsystems are locked into an active regime, while the rest present complex on-off changing dynamics, and able to simultaneously develop steady states, regular oscillations and chaotic transitions, may help to better understand cytological phenomena and to reinterpret them in a closer to reality way.

Dissipative Metabolic Networks model
DMNs consist of a set with n interconnected elements, called metabolic subsystems or MSbs. Each subsystem represents a group of enzymes aggregated in clusters and dissipatively structured. This enzymatic set is considered as an individual catalytic entity.
The subsystems receive both input fluxes (the substrates of the biochemical reactions) and regulatory signals, which may be of three types: activatory, inhibitory and total inhibitory (all-or nothing type). These interactions correspond with the activity of allosteric enzymes and regulatory enzymes of covalent modulation.
Each MSb converts the input fluxes and regulatory signals into a unique or several output activities and sends one or several fluxes and regulatory signals to other subsystems.
The output activity can be periodic or stationary (In accordance with the experimental observations, in which most of the patterns are stationary or correspond with periodic oscillatory behaviors [33]).
The conversion from input to output activity is performed in two stages. In the first one, the input fluxes are transformed in an intermediary activity using some of the ''flux integration functions''. In the second stage, the received regulatory signals originate a ''regulatory signal integration'' which varies the intermediary activity.
The magnitude of the influence of each regulatory signal is defined by an associated regulatory coefficient.
In order to simplify the assumptions, we will take first the periodic oscillations to be harmonic; these harmonic oscillations will be interspected later with transition regimes that are combinations of the harmonic oscillations preceding and following the transition. This will result in nonlinear oscillations with different levels of complexity (see, figure 6d), as could be expected in the cellular conditions ''in vivo''.
So the activity of the i-th subsystem will take the form y i (t) = (A 0 ) i +A i sin(v i t), where (A 0 ) i is the mean amplitude, A i is the amplitude of the oscillation and v i is the frequency. Since y i (t) has to be non-negative, we will take (A 0 ) i $0 and 0#A i #(A 0 ) i . We will also suppose that the values of (A 0 ) i and v i are bounded, that is, there exist constants (A 0 ) max and v max such that (A 0 ) i #(A 0 ) max and v i #v max for all i. The function y i (t) can be characterized by three variables, x i,1 , x i,2 and x i,3 , with values between 0 and 1, by taking (A 0 ) i = x i,1 (A 0 ) max , A i = x i,2 (A 0 ) i and v i = x i,3 v max (in fact, the variables x i,j may be considered as parameters, because they change very slowly in time). Hereinafter, we will identify the activity of the i-th subsystem during the kth dynamical regime with the triple x k i,1 ,x k i,2 ,x k i, 3 , and the activity of the complete set of subsystems with the array x k i,j , where i varies from 1 to n, being n the number of subsystems and where j varies between 1 and 3; the superindex k indicates the current iteration.
The activity will be steady when x i,2 = 0 or x i,3 = 0, and the subsystem will be in an inactive state when x i,1 = 0.
The normalized activity of each MSb is x k i,j . The temporal period during which each oscillation is maintained will be a constant value that we will call T. During the k-th iteration the time t will vary between (k-1)T and kT.
We proceed now to describe the way the subsystems are interconnected and how the conversion of the input activities into a unique output activity is realized. As we said above, the activity of the whole DMN in the k-th iteration will be where the index j in x k i,j define A 0 , A and v when j is 1, 2 and 3, respectively. This matrix describes the state of the whole net in each iteration and we will call it the state matrix. The integration process that transforms x k in x k+1 consists of two stages that will be fully described in which follows. The first stage is the flux integration, in which x k~xk i,j is transformed in an intermediary activity x 0k~x0

Flux integration
Each subsystem receives an input flux from a subset of the remaining ones. In the simplest case, the i-th subsystem receives flux from only one subsystem, say the l-th. The input flux will consist of the three values x k l,1 ,x k l,2 and x k l,3 , and the integrated input flux will consist of the numbers The qualitative behaviors of the amplitude and frequency in the glycolytic subsystem were studied in detail by Goldbeter and Lefever [80] (see figure 8). We have used piecewise linear functions approximating the nonlinear functions obtained by Goldbetter and Lefever. Thus, we have calculated the integrated input flux by using the functions F 1 ,F 2 ,F 3 defined as follows: p, i f xw0:9: Where p is a parameter associated to each flux integration function.
In the most general case, the i-th subsystem receives flux from more than one subsystem. We will define the integrated flux as the arithmetic mean of the received integrated fluxes.
Each subsystem can also have an outer input flux, consisting of three fixed values x 1 i,1 ,x 1 i,2 and x 1 i,3 that are integrated with the same functions F 1 ,F 2 ,F 3 that were defined before, but for other parameter values p 0 i,1 ,p 0 i,2 and p 0 i,3 . These integrated outer fluxes must be taken into account in the evaluation of the arithmetic mean together with the integrated inner fluxes.

Signal-regulatory integration
When the flux integration has been done, we proceed next with the signal regulatory integration. First, we will describe the process in the simple case in which a MSb has only one regulatory signal. Let us suppose that the i-th MSb receives a regulatory signal from the l-th one. Let x k l,j be the j-th component (which describe A 0 , A or v) of the activity of the subsystem which regulates at the beginning of the k+1-th iteration, and x 0 i,j k the j-th component of the activity of the regulated MSb after the input flux integration.
If the signal is of inhibitory type (negative allosteric modulation), then the activity x 0 i,j k is reduced by multiplying it by a factor q 0 i,j in the unit interval [0,1]. It will be 1 when x k l,j is 0, and a parameter q i,j when x k l,j is 1. In the remaining cases it will vary affinely with x k l,j , so that q 0 i,j~q i,j {1 À Á x k l,j z1 and x kz1 i,j~q 0 i,j x 0 i,j k . Let us stress that when x k l,j is 0, x 0 i,j k is unaltered, and when x k l,j is 1, the inhibition is maximum and the activity is reduced by a factor of q i,j .
If the signal is of activatory type (positive allosteric modulation), the function is similar to the one described above, once we take into account that to augment a value x in [0, 1] is equivalent to reducing 1-x to a smaller value x9 in [0, 1]. So, we will define the activatory functions by We will call the parameters q i,j ''influence coefficients''; they represent the influence of each modulator on each subsystem.
If the signal is of total inhibition type, then we will use a parameter d, that will be called threshold. The activity x 0 i,j k will be unaltered when the activity x k l,j does not reach that threshold and will be 0 when it trespass this value, that is, In the general case in which a MSb receives regulatory signals from several subsystems, they act sequentially, so that if there are r regulatory signals for each MSb, then we have r intermediary states x s,k , where s varies between 1 and r, (besides of the intermediary state x9 k obtained after the flux integration was done). In this case, we replace in the previous formulas x 0 i,j k , by x s,k i,j so that if the signal is inhibitory then we have and if the signal is activatory then we have and, finally, if the signal is of the total inhibition type (enzymatic covalent regulation) then we have ( At the end we obtain x r,k = x k+1 . Now we will describe some additional rules for the integration of the activities of each subsystem. When after the flux integration stage x 0 i,j k~0 (that is, if the input flux is 0), then we will ignore the signal regulatory integration stage and we put simply x kz1 i,j~0 . Also, if in the l-th stage of the signal regulatory integration a total inhibition acts and the threshold is reached then we will ignore the subsequent regulatory signals for that component (A 0 , A or v) and we will take x kz1 i,j~0 . Also, we will ignore regulatory signals from an inactive MSb: if the i-th subsystem receives a regulatory signal from the l-th one and x k l,1~0 , then the corresponding regulatory integration will not be performed.
Finally, when x k i,1~0 , the MSb will be off independently of x k i,2 and x k i, 3 . In this case, we will take x k i,2~x k i,3~0 . The DMNs are discrete and deterministic dynamical systems. Once the parameters of the net are fixed, there is only one solution in the net. In the simplest case, each metabolic subsystem presents a unique dynamic behavior, which can be a steady state or a periodic oscillation. However, the net can self-organize spontaneously in such a way that each metabolic subsystem shows a solution characterised by presenting uninterrupted transitions between various behaviors periodic and/or stationary. In the most complex cases deterministic chaotic solutions can emerge spontaneously, originating infinite transitions between different periodic and/or stationary behaviors. In this situation, each metabolic subsystem modifies uninterruptedly its enzymatic activity in such a way that it changes between different periodic and/or stationary regimes that never repeat its activity along arbitrarily long periods of time.
The dynamic behaviors that emerge spontaneously in the DMNs have their origin in the regulatory structure of the feedback loops, in the non-linearity of the constitutive equations of the system and in the complexity of the own dynamics that are generated in the network.

Representation of the activity of the metabolic subsystems
Now we will describe the way we pass from a sequence of terns x k 1 ,x k 2 ,x k 3 À Á representing the mean amplitude, amplitude and angular speed of a given subsystem to a continuous function representing the activity level of the subsystem. We consider a number N of transitions and, in the k-th stage, we suppose that the oscillation is harmonic, that is, the activity of the subsystem is described by a function of the form y(t) = (A 0 )+A sin(vt), where A 0~x k 1 A 0max , A~x k 2 A 0 and v~x k 3 v max and where A 0max and v max are given fixed parameters independent of the stage number and of the subsystem. The duration of the harmonic oscillation is a given parameter T h independent also of the stage and of the subsystem. In between two stages, a mixed transition regime is maintained with a duration T tr independent of the stage number and of the subsystem. If the transition goes from the k-th stage to the (k+1)-th stage then, during the T tr seconds of transition regime, the activity is given by a function of the form y(t) = A(t)y 1 (t)+B(t)y 2 (t), where y 1 (t) is the activity corresponding to the prolongation in time of the previous harmonic activity in the k-th stage, and y 2 (t) is the back-propagation in time of the subsequent harmonic activity in the (k+1)-th stage. The numbers A(t) and B(t) depend on time and indicate the weights with which the activities of the subsystem in the previous and posterior stage are present during the transition time. At the beginning of the transition, say at t = t 0 , A(t 0 ) is 1 and B(t 0 ) is 0, and at the end of the transition, say at t = t 1 , A(t 1 ) is 0 and B(t 1 ) is 1. At the rest of the transition times A(t) and B(t) vary affinely. Thus, A t ð Þ~t . Putting all this together, during the transition time the activity is given by The transition regimes are combinations of two harmonic oscillations with nonconstant coefficients A(t) and B(t) depending on time. Thus, the introduction of these transition regimes provokes the emergence of nonlinear oscillatory behaviors, both simple and complex (see for example figure 6d).

Dependence of the past
In an attempt to take into account a influence of the past in some of the DMNs, taking into account that each activity of the subsystem (besides of depending on incoming fluxes and regulatory signals) depends also on some activities developed in the past, we have considered the following aspects: the state x 00 i,j k resulting from the flux integration and signal-regulatory-integration stages is augmented or diminished in function of x k i,j and x k{1 i,j (i.e., of the state in the k-th stage and of the previous state); when In the second case, the difference 1{x 00 i,j k is reduced multiplying it by an analogous factor b9 belonging to [0, 1] whose value is 1 when x k{1 i,j {x k i,j~0 , the same constant b of the precedent discussion when x k{1 i,j {x k i,j~1 and have an affine variation with Then, To overcome the need to know two initial states in the first stage, that is, when k = 1, we have obtained the state x 2 from x 1 without taking into account the dependence of the past, which is equivalent to define x 0 = x 1 .
In some of our studies we have considered b as one of the control parameters.

Example 1
We will consider the simple DMN formed by two subsystems arranged in series with two feedback loops of regulatory signals. The MSb1 is activated by the second subsystem and the MSb2 is totally inhibited by the first subsystem when this one reaches a determinate threshold value (figure 1). The MSb1 input flux value is x As we said before, we will not consider the influence of the past for the calculation of the next state, and we will use only the flux integration and the regulatory integration. The details of these integrations will be omitted for brevity in this first stage, but the calculations will be completely developed in the next stage. The second state that we obtain is Finally, in the DMN the first MSb will fall into a single active state, corresponding to a periodic oscillation with A 0 = 0.225, A = 0.125, v = 0.35, and the second MSb is locked into an inactive state.

Example 2
Next, we have considered a DMN formed by two subsystems arranged in series, which represent in a very simplified way the main interconnections between the pyruvate dehydrogenase complex (PDH) and Krebs cycle, so-called tricarboxylic acid cycle (TCA) (Fig. 9).
Pyruvate dehydrogenase is a large multienzymatic subsystem containing many copies of three enzymes: pyruvate dehydrogenase, dihydrolipoyl transacetylase and dihydrolipoyl dehydrogenase. This metabolic subsystem (MSb1) receives an input flux of pyruvate yielding acetyl-CoA, H + and NADH by a process called pyruvate decarboxylation which is mainly inhibited by the ratio ATP/ADP [62]. When ATP, the energy-rich end product of the tricarboxylic acid cycle and oxidative phosphorylation, accumulates to high levels, the rate of formation of acetyl-CoA is slow.
Acetyl CoA, a product of the pyruvate dehydrogenase reactions is a central compound in metabolism with several functions: as input to the Krebs cycle, where the acetate moiety is further degraded to CO2, and as donor of acetate for synthesis of fatty acids, ketone bodies, and cholesterol.
One of the factors controlling the MSb2 is the NADH. The citrate synthase enzyme is the primary control point in the Krebs cycle and is negatively regulated by NADH. Isocitrate dehydro- Figure 9. Network with two metabolic subsystems: pyruvate dehydrogenase and Krebs cycle. DMN formed by two subsystems arranged in series, which represent in a very simplified way the main interconnections between the pyruvate dehydrogenase complex (PDH) and Krebs cycle. doi:10.1371/journal.pone.0003100.g009 genase and a-ketoglutarate dehidrogenase are also strongly inhibited by the negative allosteric modulator NADH [62].
Therefore, the catalytic dissipative element MSb1 is inhibited by the ATP from the second subsystem and the MSb2 receives an inhibitory signal of NADH from the MSb1.
The parameters considered in this DMN are the following ones: the first subsystem MSb1 receives an input flux of pyruvate with x 1 1,1~0 :55, x 1 1,2~0 :72 and x 1 1,3~0 :88, being the values of the parameters of the integration functions of (p 1,1 = 0.54, p 1,2 = 0.60, p 1,3 = 0.78) and b = 0.2. The initial states are the same that the ones belonging to the DMN shown in the example 1.
In order to simplify we have equalized the three influence coefficients of A 0 , A and v belonging to each inhibitory signal (q i1 = q i2 = q i3 = q i ).
First, we have considered as control parameter q 1 , the ATP influence coefficient, and we have fixed q 2 = 0.01 which represents the influence coefficient value of NADH (small values of the influence coefficient represents high inhibitory activity so for q 1 = 0.1 and q 2 = 0.01 both subsystems are strongly inhibited).
After the numeric integration of the DMN it can be observed that for 0.1#q 1 #0.9 both subsystems are active and the two subsystems exhibit one periodical activity. Some of the values of the amplitude and oscillatory frequencies are shown below: When the inhibition provoked by the ATP descends (by increasing the values of the influence coefficient) the activity of the MSb1 increases as well as the contribution of NADH and as a consequence of it the activity of the MSb2 descends.
We have considered now p 2 as control parameter, the NADH influence coefficient, and we have fixed the ATP influence coefficient value (q 1 = 0.01). This q 1 value represents a strong inhibition of the pyruvate dehydrogenase complex due to high levels of the energy-rich end product of the tricarboxylic acid cycle and oxidative phosphorylation. Under these parametric conditions the metabolic network can present different dynamic behaviors for q 2 values.
When q 2 ranges from 0.1 to 0.87, both subsystems are active and exhibit a single periodical activity. For example: As it can be observed, for q 2 = 0.1 the two subsystems are strongly inhibited and both show a low activity. When increasing p 2 values (0.1,d#0.87) the inhibition provoked by the NADH decreases and as a consequence of it the activity of the MSb2 is increased and the MSb1 activity diminishes due to the increment in the production of ATP.
For 0.87,q 2 #1 the two subsystems undergoes a variation of the oscillatory evolution, emerging very complex transitions and the output activity of the both subsystems make uninterrupted transitions between more of 200 different kinds of periodic oscillations.
Under experimental conditions, the complex temporal structure in Krebs cycle has been observed [78] and it has been verified the complex transitions that characterize the intermediates of the Krebs cycle. See figure 5a for more details.
In the network the two metabolic subsystems remain always active for all the parametric conditions studied.

Experimental calcium oscillations in Xenopus laevis oocyte
In order to study chaotic oscillations in a subsystem with an onoff changeable regime we carried out experimental observations of intracellular calcium oscillations in Xenopus laevis oocyte.
Membrane currents were recorded with a standard twoelectrode voltage clamp (Warner Instruments, Oocyte Clamp OC-725C) and digitized in a PC computer (Digidata 1200 and Axoscope 8.0 software, Axon Instruments). Oocytes were continually superfused with Ringer's solution at room temperature (,22uC).
The membrane was usually voltage clamped at 260 mV. Fetal Bovine Serum (FBS) (Sigma-Aldrich) diluted 1:1000 in Ringer solution were perfused to achieve the generation of voltage oscillations.