INTRODUCTION

TOP ABSTRACT INTRODUCTION THE MODELS RESULTS DISCUSSION REFERENCES

Under some conditions, the concentrations of metabolites of the glycolytic pathway oscillate (Duysens and Amesz, 1957; Betz and Chance, 1965; Hess and Boiteux, 1980). In suspensions of intact yeast cells that have been harvested during the growth-phase transition between using glucose and ethanol as the substrate, then starved and subsequently presented with glucose and cyanide, these oscillations are sustained (Richard et al., 1993, 1996a). The reproducible observation of limit-cycle oscillations in well-defined, intact cells, has made it possible to study intracellular dynamics and its control under stationary conditions (Hess and Boiteux, 1980; Termonia and Ross, 1981; Richard et al., 1996a,b; Reijenga et al., 1998).

For the oscillations to be observable, they had to be synchronous in the billions of cells in the cuvette. Indeed, mixing of cell populations that oscillated 180° out of phase, only caused a transient disappearance of the oscillations (Richard et al., 1996b), evidencing some sort of active synchronization. Acetaldehyde synchronizes the oscillations between the cells (Richard et al., 1994, 1996b). This identification has made this system also an attractive model for the study of the control of intercellular dynamics (Wolf and Heinrich, 2000; Bier et al., 2000).

For the mechanism of these oscillations, a number of proposals exist. Some focus on phosphofructokinase as the pacemaker, others on the stoichiometry of the reaction network (Sel'kov, 1975; Goldbeter, 1996). The notion that control of the oscillations would have to be confined to phosphofructokinase has been shown to be unjustified for either type of model (Teusink et al., 1996; Bier et al., 1996). Preliminary results suggest that control of the frequency of the oscillations is distributed over a number of molecular processes, one of which includes glucose transport (Reijenga et al., 1998; cf., Markus et al., 1984).

A crucial issue in how a dynamic cellular system controls itself is the active networking in which its macromolecules engage. Upon a perturbation anywhere in the system, the enzymes around that perturbation respond in terms of a change in rate. The changes in rates result in secondary perturbations, which again result in rate changes of surrounding enzymes. The key property of stable systems is that the resulting parallel chains of changes precisely return the system to its original state (Westerhoff and Van Dam, 1987). To understand the internal regulation of cellular systems, it is important to understand how the various changes match precisely (Kahn and Westerhoff, 1993). In systems at steady state, such an analysis is difficult. Because cause and effect cannot be distinguished temporally; one should analyze transient perturbations, or cut up regulatory links in the system (Snoep et al., 1999). For systems engaged in limit cycles, it is more feasible to trace the paths of the system dynamics, because these allow the analysis of phase (Betz and Chance, 1965) and amplitude (Richard et al., 1996a) differences. Because of the possibility to measure reliably and rapidly the glycolytic intermediates during the oscillations (Richard et al., 1996a), we set out to determine the chain of events in synchronizing yeast glycolytic oscillations. Such a study should not only address the chain of events within individual cells, but also that between the cells.

Surprisingly in this respect, Richard et al. (1996a) found oscillations in glycolytic metabolites to be virtually confined to the hexose phosphates in the upper part of the glycolytic pathway. Oscillations were observed in the concentrations of glucose 6-phosphate (G6P), fructose 6-phosphate (F6P), fructose-1,6-bisphosphate (F16P). The concentrations of dihydroxyacetone phosphate (DHAP), glyceraldehyde 3-phosphate (GAP), 3-phosphoglycerate (3PG), phosphoenolpyruvate (PEP) and pyruvate oscillated at a relative amplitude of less than 2% for DHAP and less than 30% for the other metabolites. The concentrations of 1,3-biphosphoglycerate (1,3BPG) and 2-phosphoglycerate (2PG) were below the detection limit. These findings suggested that neither the oscillations in intercellular acetaldehyde concentration, nor the cell-cell synchronization could be understood in terms of a cause-effect chain running through the intermediates of the glycolytic pathway itself.

To rationalize this, Richard et al. (1996a) pointed out that the coenzymes engaged in glycolysis, i.e., ATP(ADP) and NAD(H), also oscillated in terms of their concentrations. Following this lead, we here elaborate a mechanism for the propagation of the oscillations and their phase through the glycolytic network in and between the individual cells. By suppressing other possible modes of transfer of dynamics, we show that this mechanism may indeed be responsible for the experimental observations in dynamically coupled oscillating yeast cells.

	THE MODELS

TOP ABSTRACT INTRODUCTION THE MODELS RESULTS DISCUSSION REFERENCES

We have set up a minimum model of glycolysis in yeast that suffices to describe the essence of the experimental observations of Richard et al. (1996a,b). This model contains lumped reactions of the glycolytic pathway and includes production of glycerol, fermentation to ethanol and exchange of acetaldehyde between the cells, and trapping of acetaldehyde by cyanide. The system was analyzed with the software packages AUTO (Doedel, 1981) as well as a Runge-Kutta-Merson integration routine. In the situation under study, respiration is absent, and ethanol is the main end product of glycolysis. The reactions and fluxes under consideration are shown in Fig. 1. For the extracellular fluxes, we take into account a glucose flux into the cell. Ethanol and glycerol concentrations are considered fixed due to the large extracellular reservoir with which they are supposed to equilibrate. Minor side fluxes are fluxes resulting from the exchange of acetaldehyde between the cell and the surrounding medium, a glycerol production flux, and an extracellular sink for acetaldehyde. The cell membrane is treated as impermeable for metabolites other than glucose, acetaldehyde, and ethanol.

View larger version (56K): [in this window] [in a new window]   FIGURE 1   Schematic representation of the anarobic glycolytic pathway. Reactions: HK-PFK, lumped reaction of hexokinase and phosphofructokinase; ALD, aldolase reaction; GLY, glycerol-producing branch; GAPDH, reaction of glyceraldehyde-3-phosphate dehydrogenase; PGK, reaction catalyzed by phosphoglycerate kinase; PK, lumped reactions of phosphoglycerate mutase, enolase, and pyruvatekinase; PDC, reaction catalyzed by pyruvate decarboxylase; ADH, alcohol dehydrogenase reaction; ATPase, total cellular ATP consumption; CYA, degradation of acetaldehyde by cyanide. Fluxes: glucose influx, exchange of acetaldehyde between the cell and the external medium.

The intracellular reactions that are taken into consideration are: ALD, aldolase; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; PGK, phosphoglycerate kinase; PDC, pyruvate decarboxylase; ADH, alcohol dehydrogenase. The model includes the degradation of the coupling substance in the external solution (CYA), i.e., the reaction of acetaldehyde with cyanide. Consumption of ATP is taken into account by including an ATPase. Some reactions describe lumped processes. This concerns the reactions catalyzed by hexokinase (HK) and phosphofructokinase (PFK) (lumped in reaction HK-PFK), the glycerol-producing steps (resulting in reaction GLY) as well as phosphoglycerate mutase, enolase, and pyruvate kinase (lumped together in reaction PK).

With the exceptions of GAPDH and PGK, all reactions are modeled as if irreversible. The reactions catalyzed by GAPDH and PGK are near their equilibria with the equilibrium constants q_GAPDH = k_GAPDH+/k_GAPDH = 0.0056 (Byers, 1982) and q_PGK = k_PGK+/k_PGK = 3225 (Bergmeyer, 1974). The balance equation for 1,3-bisphosphoglycerate (1,3BPG) reads

(1a)

with

(1b)

and

(1c)

[PTP], [1,3BPG], and [3PG] denote the concentrations of the pools of triosephosphates (dihydroxyacetone phosphate and glyceraldehyde-3-phosphate), 1,3BPG and 3PG, respectively. The equilibrium constants for GAPDH and PGK lead to a very small steady-state concentration of 3PG and, hence, to a short average life time as compared to the other metabolites. This warrants a quasi-steady-state approximation for [1,3BPG] (see Heinrich and Schuster, 1996). From d[1,3BPG]/dt = 0, it follows that v_GAPDH = v_PGK, and, for the concentration of 1,3BPG

(2)

In this way, the number of variables was reduced by one. The rate equation for the GAPDH-PGK reaction reads

(3)

As a consequence, the GAPDH-PGK reaction can be understood as quasi-trimolecular. The resulting model is the 9-variable system shown in Fig. 2. The variables consist of the concentrations of the following metabolites: S₁, glucose; S₂, fructose-1,6-bisphosphate; S₃, pool of the triosephosphates, glyceraldehyde-3-phosphate, and dihydroxyacetone phosphate; S₄, 3-phosphoglycerate; S₅, pyruvate; S₆, acetaldehyde in the cell; S₆^ex, extracellular acetaldehyde; A₂, ADP; A₃, ATP; N₁, NAD⁺; and N₂, NADH. In the following, we consider the cellular pools of the adenine nucleotides ADP and ATP as well as that of the nicotinamide adenine dinucleotides NAD⁺ and NADH as conserved moieties with

(4a)

and

(4b)

For the mathematical description, simple rate laws are used for all enzymes. We describe the activities of the enzymes with linear or bilinear terms. Only for the HK-PFK reaction, regulatory properties are taken into account. PFK is a strongly regulated enzyme and has several effectors. It is activated by F6P, F16P, ADP, and AMP, whereas ATP leads to inhibition. Here, the inhibition by its substrate ATP is considered by a factor
where K_i and n denote the inhibition constant and the cooperativity coefficient of that regulation (Heinrich and Rapoport, 1975).

View larger version (50K): [in this window] [in a new window]   FIGURE 2   Schematic representation of the models. The 9-variable model includes all shown metabolites as variables, whereas, in the 6-variable model S3, S4, and S5 are treated as parameters. Variables: S1, glucose; S2, fructose-1,6-bisphosphate, S3, pool of the triosephosphates, glyceraldehyde-3-phosphate, and dihydroxyacetone phosphate; S4, 3-phosphoglycerate; S5, pyruvate; S6, acetaldehyde in the cell; S6ex, extracellular acetaldehyde; A3, ATP, and N2, NADH.

The glucose influx into the cell is assumed to be constant, whereas the flux of S₆ across the plasma membrane is taken to be proportional to its concentration difference across that membrane

(5)

(6)

where  acts as a coupling parameter (A, cell surface area; V_c, cellular volume; P, permeability). The ratio of the cellular and the extracellular volume (V_ex) is denoted by . Assuming a homogeneous distribution of the metabolites in the intracellular and in the extracellular solution, the differential equation system of the model follows:

(7a)

(7b)

(7c)

(7d)

(7e)

(7f)

(7g)

(7h)

(7i)

with the rate equations

(8a)

(8b)

(8c)

(8d)

(8e)

(8f)

(8g)

(8h)

(8i)

One rewrites Eq. 3 for v_GAPDH-PGK in terms of the variables

(8j)

	RESULTS

TOP ABSTRACT INTRODUCTION THE MODELS RESULTS DISCUSSION REFERENCES

Description of the experimental observations: the 9-variable model

We first examined whether the 9-variable system (Eqs. 7a-i and 8a-j) could describe the dynamics we observed experimentally (Richard et al., 1996a,b). The reference parameter set shown in Table 1 was chosen, such that the steady-state concentrations and fluxes were in an experimentally realistic range. In the following, we use the degradation rate constant of acetaldehyde (k₉) as the bifurcation parameter. Experimentally, k₉ can be modulated by changing the cyanide concentration. For each value of k₉, the system had one steady state with nonnegative metabolite concentrations. This steady state was stable if the degradation rate k₉ was small and unstable for higher values of this parameter (see Fig. 3). The instable steady state arose via a subcritical (http://mrb.niddk.nih.gov/glossary/glossary.html) Hopf bifurcation in the point, indicated by k_9(H) in Fig. 3. In the region k_9(S)  k₉  k_9(H), a stable steady state and an unstable limit cycle coexisted. k_9(S) denotes a saddle-node bifurcation of the limit cycle, with the consequence that, for k_9(S)  k₉  k_9(H), a stable steady state and an instable limit cycle coexisted with a stable limit cycle. In Fig. 3, the steady-state concentration and the minima and maxima of the oscillations are drawn for the variable N₂.

View larger version (14K): [in this window] [in a new window]   FIGURE 3   Bifurcation diagram for the variable N2 (NADH) in mM. The bifurcation parameter k9 is the pseudo-first-order rate constant of the extracellular degradation of acetaldehyde. Solid line, stable steady state; dashed line, unstable steady state. Empty circles, maxima and minima of unstable limit cycles; full circles, maxima and minima of stable limit cycles. Points: S, point of a saddle-node bifurcation of the limit cycle; H, point of Hopf bifurcation. Parameters as given in Table 1.

For the set of reference parameters and k₉ = 80.0 min¹, we analyzed the oscillations in greater detail. The oscillation period was T = 0.141 min, i.e., approximately six times shorter than the experimental one (Richard et al., 1996a). Decreasing all rate constants brought the period closer to the experimental one, at the cost of the model average flux becoming lower than the experimental flux. Introducing saturability of the rate equations should increase the period while keeping the flux constant, but we preferred to keep the model simpler, not to get bogged down in overparametrization. After all, the focus of this paper is to demonstrate the feasibility of a mechanism of the propagation of dynamics, rather than an exact fitting of experimental data.

Figure 4 shows an integration in time of NADH (N₂) and ATP (A₃) concentrations. The two metabolites had a phase shift (152°), which was comparable to that in the experiments. Also, the phase shift of S₆^ex and N₂ of 138° was in accordance with the experimental situation (Richard et al., 1996a). The phase relationship between S₂ and N₂ was not quite reproduced by this core model. In Table 2, the mean values and relative amplitudes of all variables are given. The mean values correspond to the experimental results shown in Richard et al. (1996a). The concentrations for internal and external acetaldehyde are also in agreement with a paper of Stanley and Pamment (1992), where an accumulation of acetaldehyde inside of cells was demonstrated.

View larger version (15K): [in this window] [in a new window]   FIGURE 4   Integration of NADH and ATP of the 9-variable model in time. The period of the oscillation is T = 0.141 min, phaseshift of NADH (thin line) and ATP (thick line) equals 152°. Parameters as in Table 1, and k9 = 80.0 min1.

We did not find the same strong separation in the relative amplitudes of the oscillating metabolites as was observed in the experiments. This may again be due to the fact that the enzyme kinetics in the model are simpler than those in reality. For the parameter set we used, the relative amplitudes in S₁ and A₃ were very high, that in S₅ was very small, and, for all other metabolites, the relative amplitude was moderate. The essence of an observation made by Richard et al. (1996a) was reproduced: Following their rule of thumb that the relative amplitude (a) of a driven oscillation should be smaller than the amplitude of the driving oscillation, the oscillations in S₆ cannot be explained in terms of the oscillations being transduced only through the backbone of the glycolytic pathway. For, a_S4 a_S3 and a_S6 > a_S5. Richard et al. noted that the oscillations in the ATP/ADP ratio did exceed the oscillations in S₄ in terms of relative amplitude. The former also exceeded those in the NADH/NAD ratio, which, in turn, exceeded those in S₆. This is also observed in our model calculations.

Feasibility of the proposed intracellular propagation mechanism; reduction to the 6-variable model

A numerical model has the disadvantage that it may not precisely mimic reality, but the advantage is that it is free of experimental error. Consequently, the finding that, also in a numerical model, substantial oscillations in acetaldehyde arose in the absence of such oscillations all along the backbone of glycolysis suggested to us that this observation was not due to experimental error, hence worthy of a nontrivial explanation. In coming up with the explanation, we were inspired by the rule of thumb that the relative amplitude of a driven oscillation should be smaller than that of the oscillation driving it (Richard et al., 1996a). We propose that the oscillations originating around phosphofructokinase engages the adenine nucleotides and especially their concentration ratio (ATP/ADP) in an oscillation of substantial amplitude. Coupling between the oscillations in ATP hydrolysis free energy and the oscillations in the NADH redox free energy is proposed to be provided by the GAPDH and PGK reactions. Further down glycolysis, the alcohol dehydrogenase reaction should likewise couple oscillations in the NADH/NAD ratio to oscillations in the concentration of acetaldehyde (cf. Richard et al., 1996a). The oscillations in acetaldehyde should then mediate intercellular synchronization, implying that, in the neighboring cells, there is an influence from the bottom of glycolysis (at the alcohol dehydrogenase reaction) back up, ultimately to phosphofructokinase. The validity of the rule of thumb can only be proven for systems in the Onsager realm (Richard et al., 1996a), which the present system is not. Consequently, the feasibility of this mechanism for the propagation of the dynamics within the cell requires proof. We first set out to test the hypothesis that the oscillating signal can be transduced through the pathway, through coupling of the oscillations in NADH redox free energy and oscillations in the Gibbs free energy of ATP hydrolysis, bypassing the carbon metabolites of the lower part of glycolysis. To this aim, the concentrations of the carbon metabolites in the lower part of glycolysis were fixed in the model. If the above hypothesis is correct, this should not interfere with the oscillations.

The differential equation system for the resulting 6-variable model reads

(9a)

(9b)

(9c)

(9d)

(9e)

(9f)

with Eq. 8a-j for the rates, and treating S₃, S₄, and S₅ as (fixed) parameters. By fixing S₃, S₄, and S₅ to the steady-state concentrations ₃(k₉), ₄(k₉), and ₅(k₉) calculated from the 9-variable model, it was ensured that the pathway remained balanced as a whole in the stationary state. Accordingly, the stationary concentrations of the variables in the reduced model were identical to those in the 9-variable model.

The stability of the steady state was analyzed by calculating the eigenvalues of the Jacobian matrix of the 6-variable model. For the set of reference parameters, the result is shown in Fig. 5 for N₂. Again, the steady state was stable for small degradation rates and unstable for higher ones. The region of stability exceeded that of the 9-variable model (see Fig. 5).

View larger version (11K): [in this window] [in a new window]   FIGURE 5   Steady-state concentration of NADH in the 6-variable model as a function of k9. Thick line, stable steady state; thin line, unstable steady state. Parameters as in Table 1.

For the degradation rate k₉ = 80.0 min¹, we investigated the dynamics of the reduced model and compared it to the corresponding behavior of the 9-variable model. By integration, we confirmed that the stationary state was unstable for that parameter set. Moreover, we found a stable limit cycle with an oscillation period of T = 0.104 min. In Table 3, the corresponding mean values and relative amplitudes of the oscillating metabolites are given. The mean values were in the same order of magnitude as in the 9-variable model. The relative amplitudes were somewhat reduced compared to the ones in that model, but were still much higher for S₁ and A₃ than for the other metabolites. The phase shifts of NADH and ATP amounted to 146°, that of NADH and external acetaldehyde to 162°. The phase shifts were not affected much by the model reduction.

Importantly, S₆ oscillated due to the primary oscillations around reaction v₁: oscillatory dynamics had traveled from v₁ to S₆ at fixed concentrations of S₃, S₄, and S₅. This proves that the dynamics need not spread through the backbone of glycolysis, but can do so through the ATP/ADP and NADH/NAD coenzymes. In addition, the rule of thumb of relative amplitudes continued to apply (cf. Table 3).

The fixation of S₃, S₄, and S₅ to their steady-state values caused an imbalance of the time-averaged fluxes beyond the steady state, i.e., in the case of oscillations. We therefore also examined a case in which S₃, S₄, and S₅ were fixed such that the time averaged fluxes were balanced. This required

(10a)

(10b)

and

(10c)

where T is the oscillation period. When these conditions 10a-c were fulfilled, the oscillations persisted. For the set S₃ = 0.53393 mM, S₄ = 0.56646 mM, S₅ = 8.35040 mM, for instance, the modulus of all above integrals was smaller than 0.05% of the average flux through the pathway.

The transduction of the dynamics between the cells

We also asked if, in the absence of oscillations in the glycolytic backbone, the cells would be able to synchronize. We considered two cells that were taken identical with respect to the stoichiometry, kinetics, and permeability properties. Two cells at different phases were made to share their extracellular acetaldehyde concentration. The time integration is shown in Fig. 6. The initial phase shift vanished and the cells synchronized, demonstrating that synchronization could be accounted for by the model. For further investigations of the intracellular coupling of glycolytic oscillations in yeast cells, see Bier, et al. (2000) and Wolf and Heinrich (2000).

View larger version (20K): [in this window] [in a new window]   FIGURE 6   Integration of two interacting cells. Shown are the NADH concentrations in both cells. The initial phase shift vanishes and the cells synchronize in time. Parameters as in Table 1, and k9 = 80 min1, and S3, S4, S5 fixed to the steady state values. A, B, C, and D are subsequent snapshots, i.e., after 0, 1.5, 3.5, and 29.5 min, respectively.

	DISCUSSION

TOP ABSTRACT INTRODUCTION THE MODELS RESULTS DISCUSSION REFERENCES

We have compared oscillations in a model for glycolysis to a model in which the metabolites in the lower part of glycolysis are fixed. Both models are core models, i.e., they have a limited number of variables and are intended to study essential features of the glycolytic oscillations. The parameters were adapted to reproduce experimentally observed average fluxes and concentrations. The amplitudes and phases of the oscillating metabolites matched the experimental data qualitatively. The period of the oscillations was shorter than was observed experimentally, but this may be related to the absence of saturability in most rate equations.

Using the second model, we have demonstrated that oscillations in acetaldehyde can occur in a model where the concentrations of the carbon metabolites in the lower part of glycolysis are fixed. This result shows that it is possible to transduce the oscillations from the oscillation nucleus of the model (i.e., the HK-PFK reaction) to acetaldehyde solely via coenzyme oscillations. The latter consisted of oscillations in the ATP/ADP ratio and oscillations in the NADH/NAD ratio. The coenzyme oscillations were coupled through GAPDH and PGK. By engaging alcohol dehydrogenase, the oscillation in NADH/NAD led to an oscillation in extracellular acetaldehyde concentrationsbypassing the lower part of the carbon skeleton of the pathway.

For biology, it is important that intracellular reaction networks engage in functional behavior. Often this means that a steady flux or concentration should be produced. In other cases, a steady oscillation is functional. In almost all cases, the dynamic behavior should be robust, i.e., not be abolished by spurious challenges deriving from fluctuations in external parameters or internal variables. Such robustness is achieved by virtue of the sensitivity (elasticity (Burns et al., 1985)) of the processes to the concentrations of the metabolites. The networking of the intracellular concentrations via the molecular processes returns the system to its original state after a temporary perturbation (Westerhoff and Van Dam, 1987), or to a nearby state upon a small sustained modulation. The networking cannot be readily determined for systems at steady state, but a system with oscillations offers unique possibilities, both in terms of relative amplitudes and in terms of phase relationships. For the system of yeast glycolytic oscillations, we have here achieved an understanding of the dynamic networking. We have shown that, in essence, oscillations arising around phosphofructokinase propagate through the system via the ATP/ADP ratio, the NADH/NAD ratio and the intracellular acetaldehyde concentration to the intercellular acetaldehyde concentration. In addition, we have demonstrated that the dynamics propagates into other cells in the same suspension to the extent of synchronizing them. Of course, the demonstration was for the mathematical replicon of the experimental system only, but the similarities between model and experiments suggest that the proposed mechanism of propagation of the dynamics also applies to the yeast cell suspension. Thus, we may have achieved the first demonstration of the traveling of dynamics down a chain of intra- and intercellular processes.

Our approach of fixing the concentrations of some metabolites may be used in general to distinguish metabolites that are important for an oscillation (masters, or communicators) from metabolites that fluctuate as the result of an oscillation that originates elsewhere (mere slaves). Fixing the latter should have no effect on the oscillation itself. The concentrations must be fixed in such a way that the average fluxes remain balanced during the oscillation. It is not unconceivable that application of the approach to models of synchronizing intracellular calcium oscillations may also enhance the understanding of the networking in that system (e.g., Höfer, 1999).

In an earlier investigation, it was shown that a sufficient interaction between living cells may lead to coordinated behavior. The latter may consist of synchronous or regular asynchronous oscillations (Wolf and Heinrich, 1997). For very weak or no interaction, a conservation of the initial phase shift of the two cells is expected. The earlier work was for a highly simple model that was not directly comparable to the oscillations that occur in yeast. The 9-parameter model used here is highly representative of the glycolytic pathway in yeast cells and confirmed the conclusions of the earlier work in that the coupling of two identical cells can lead to synchronization of the oscillations. The 6-parameter model then showed that the intercellular communication is possible despite the fixing of the metabolite concentrations of the lower part of glycolysis.

For high trapping rates, the cells communicated only slowly via the exchange of acetaldehyde. The individual cells continued to oscillate, but synchronization activity ceased to exist. Macroscopically, this could appear as a steady-state situation. The modeling results offer a possible explanation for the experimental finding that, macroscopically, no sustained glycolytic oscillations can be detected when the cyanide concentration exceeds a certain value in the medium (Richard et al., 1994). Indeed, our results confirm that the trapping rate of acetaldehyde plays a crucial role in the occurrence of the sustained oscillations in glycolysis. The trapping is important not only for intercellular signal transduction, but also for oscillations in a single cell. Sustained oscillations were only found for values of k₉ (acetaldehyde trapping rate) higher than a critical value, whereas, for lower values of k₉, the steady state was always stable. This corresponds to the experimental findings that macroscopic oscillations in populations of intact yeast cells can be induced by certain trapping rates of the coupling substance acetaldehyde by addition of certain concentrations of cyanide to the medium (Richard et al., 1994).

The early days of the analysis of the control and regulation of cellular processes were dominated by the oversimplification that they should be controlled by single rate-limiting steps. For metabolic fluxes at steady state, it has since become clear that biology is more sophisticated than this; control of flux tends to be distributed (Jensen et al., 1995; Fell, 1997). Although the important role of phosphofructokinase in glycolysis may have suggested otherwise, it has also been shown that this enzyme need not control the oscillations. Also, here, the situation turned out to be potentially complex in that frequency and amplitude may well be controlled to different extents by different molecular processes; there is no such thing as the controller of glycolytic oscillations (Bier et al., 1996; Teusink et al., 1996). The present paper elucidates another tier of the richness of biochemical phenomena, i.e., the fact that various steps may contribute to the propagation of oscillatory dynamics through a system and to the synchronization of the oscillations. In this case, the important steps were a bit unexpected, because glycolytic oscillations might be expected to be confined to the carbon pathway of glycolysis itself.

Control of the oscillations does not only address frequency and amplitude: the observed synchronization implies that, also, the relative phase of the oscillations of two cells became a controlled property. That biology is this complex may be unattractive from the point of view of elementary physics, but it is quite in line with the robustness required of living organisms. Also, may attract those interested in thorough and functional biocomplexity. Much experimental and theoretical work may lie ahead here (Westerhoff et al., 1999).

That the coenzymes are engaged in the oscillations is even more surprising in view of their involvement in many other aspects of cell function. The oscillation of the ATP/ADP ratio is likely to be felt by a great variety of processes. Although both ATP hydrolysis free energy and NADH redox free energy oscillate (180° out of phase), the cells did not, however, run out of free energy. During the oscillations, the energy charge remained above 0.5. The ATP/ADP ratio remained well above 0.7. Apparently, sustained oscillations in Gibbs energy are compatible with continued cell function. Termonia and Ross (1981) have predicted this from modeling results and have even suggested that oscillations may benefit the thermodynamic efficiency of the cell operations.