Abstract

The response of oscillatory systems to external perturbations is crucial for emergent properties such as synchronisation and phase locking and can be quantified in a phase response curve (PRC). In individual, oscillating yeast cells, we characterised experimentally the phase response of glycolytic oscillations for external acetaldehyde pulses and followed the transduction of the perturbation through the system. Subsequently, we analysed the control of the relevant system components in a detailed mechanistic model. The observed responses are interpreted in terms of the functional coupling and regulation in the reaction network. We find that our model quantitatively predicts the phase-dependent phase shift observed in the experimental data. The phase shift is in agreement with an adaptation leading to synchronisation with an external signal. Our model analysis establishes that phosphofructokinase plays a key role in the phase shift dynamics as shown in the PRC and adaptation time to external perturbations. Specific mechanism-based interventions, made possible through such analyses of detailed models, can improve upon standard trial and error methods, e.g. melatonin supplementation to overcome jet-lag, which are error-prone, specifically, since the effects are phase dependent and dose dependent. The models by Gustavsson and Goldbeter discussed in the text can be obtained from the JWS Online simulation database: (https://jjj.bio.vu.nl/models/gustavsson5 and https://jjj.bio.vu.nl/models/goldbeter1)

Introduction

Phase shifts induced by perturbation of an oscillatory system can be characterised in a so-called phase response curve (PRC). This graphical analysis has been applied to study the dynamics of biological oscillators in individual organisms and in populations [16], but it can also be used to quantify coupling strengths and to predict complex dynamics resulting from perturbations [7]. In behavioural studies in circadian systems [1,8], these PRCs fulfil a largely descriptive role, and although they can be used to manipulate phase shifts, for instance, to reduce jet-lag symptoms via melatonin addition [9], this is mostly a trial and error empirical approach. Although some advances towards more detailed mathematical models of the mammalian circadian clock have been made [10], melatonin as an effector is not included. With increasing knowledge of the system, more detailed mathematical models can be constructed that allow for a deeper understanding of the mechanism underlying the PRC. Typically, such studies are at the level of the individual's physiology, for instance, in cardio-respiratory studies [11], or at the cellular level in neuroscience.

For a mechanistic interpretation of the perturbation experiments, mathematical models for the system are crucial (e.g. see [6] for applications in neuroscience to study oscillations, phase shifts, and synchronisation). Another oscillatory system, for which detailed mechanistic models exist and many aspects have been studied in detail, is yeast glycolysis: e.g. the conditions necessary for oscillations and the features of oscillations in populations of cells and cell extracts [12]; oscillations in single cells [13]; acetaldehyde (ACA) acting as the communicating metabolite leading to the phenomenon of synchronised oscillating populations [14,15]; and mechanisms in the glycolytic pathway that are required for oscillations to occur (auto-catalysis [16] or the allosteric regulation of phosphofructokinase (PFK) [17,25]). In previous studies, we have also constructed and validated a detailed enzyme mechanistic model of glycolytic oscillations under a variety of conditions [13,1820]. Possibly because these studies were for a long time limited to synchronised cultures, PRC analyses have not been widely used for glycolytic oscillations in intact cells.

Experimental investigations have been conducted into the phase response of glycolytic oscillations in yeast cell extracts, to the addition of adenosine phosphates [2123] and ACA [22,24]. In [2528], it was demonstrated that the instability leading to oscillations in the glycolytic pathway and many of the features of these oscillations (including the shape of the PRC resulting from adenosine phosphate pulses) could be understood when considering a model of allosteric regulation of PFK by the adenosine phosphates inducing a switch between active (R) and inactive (T) states of the enzyme. Further stressing the role of adenosine phosphates, a model analysis of the stationary properties of oscillations, and comparison to experimental data from cell extracts and intact cell suspensions, showed that the stoichiometry of the ATP–ADP–AMP system and the allosteric regulation of PFK could be responsible for the instability leading to oscillations [29].

Since it has been shown that single cells in isolation do oscillate [13], the emergence of synchronised oscillations at the level of a population of cells must be due to a dynamic adaptation of oscillations in individual cells to an external entraining signal. It is generally accepted that the metabolite responsible for the communication between such oscillating cells is ACA [14,15]. In [14], this was investigated by mixing two out-of-phase populations and observing subsequent synchronisation and external ACA oscillations. In addition, the authors demonstrated the phase shift response of an oscillating population to external ACA addition. It was also noted in [30] that the phase shifts induced by the addition of ACA to cell suspensions depend on the phase of oscillations at the time of addition. In [15], it was shown that it is possible to entrain oscillating cell populations using oscillating external ACA concentrations and the recovery of phase and amplitude upon a 180° shift in the entraining signal was characterised. In [29], it is proposed that, due to its effect on the frequency, the NAD/NADH feedback system plays an important role in the cell-synchronisation process. However, the specific mechanism by which cells are able to dynamically adapt their oscillations to a forcing signal, whether it be from a surrounding population of cells leading to synchronisation or from an external source leading to entrainment, has not been established, and it is the subject of this study.

In [31], entrainment of oscillating isolated cells by an oscillating external cyanide concentration was shown. Depending on the phase of the cells at the time of perturbation, varying degrees of phase retardation were observed in the PRC, but no mechanistic interpretation was made. In the current study, we first used our detailed kinetic model of yeast glycolysis [20] to analyse the cyanide perturbations published in [31]. Although cyanide does trap ACA via lactonitrile formation [32], it binds also to many other metabolites [33]. We therefore extended the experiments to include direct pulses of ACA. The comprehensive experimental data set showed previously undetected phase advancement for perturbations around a specific intrinsic phase of oscillation. Our model simulations for either cyanide or ACA perturbations agreed well with the data, including the prediction of this phase advancement, and the greater likelihood to observe these in the new experiments using direct ACA pulses. A detailed modelling analysis was conducted to unravel the mechanism via which the ACA perturbation is transduced through the reaction network, eventually leading to a specific phase advancement or retardation. A key role for the PFK was revealed, controlling both the phase shift and adaptation time to an external forcing signal.

Methods

Cell preparation

Yeast cells (Saccharomyces cerevisiae, X2180) were grown in a medium containing 10 g/l glucose, 6.7 g/l yeast nitrogen base and 100 mM potassium phthalate (pH 5.0) on a rotary shaker at 30°C, and as described in [34]. The cells were harvested at the diauxic shift, centrifuged (3500g for 5 min), and washed twice in 100 mM potassium phosphate buffer (pH 6.8). Subsequently, cells were resuspended in the glucose-free potassium phosphate buffer and starved on a rotary shaker for 3 h at 30°C. Cells were then washed in a buffer solution and stored at 4°C until use.

Experimental procedure

The experimental setup and procedure were described in detail in [13]. The cells were positioned at the bottom of a microfluidic flow chamber with four inlet channels and with a cell–cell distance of 10 µm using optical tweezers. By positioning the cells far apart, cell–cell interactions, by e.g. ACA, were avoided. The flows in the channels were laminar and solutions from different channels could only mix by diffusion. The extracellular environment in the cell array area could thus be controlled both spatially and temporally by adjusting the flow rates in the different inlet channels, causing the cell array area to be covered by the chemicals from the intended channel. All channels contained 100 mM potassium phosphate buffer (pH 6.8). Glycolytic oscillations were induced by a 20 mM glucose and 5 mM KCN solution for 10 min before the perturbations. For perturbations with ACA, a solution containing 20 mM glucose, 5 mM KCN and 1 mM ACA was introduced for 20 s. The cell responses were studied by imaging the NADH autofluorescence from the individual cells using an epi-fluorescence microscope (DMI6000B, Leica Microsystems) and an EM-CCD camera (C9100-12, Hamamatsu Photonics), where images were acquired every other second.

Model simulation of experimental procedure

The Gustavsson model was previously constructed for oscillating single cells in the microfluidic environment [20]. In this experiment, the solution containing ACA contains the same cyanide concentration as the ACA-free solution (5 mM). Cyanide reacts with ACA in a lactonitrile formation process [32] leading to a decrease in the free ACA. Using the rate constant (0.035 mM−1 min−1) for this process in PBS determined in [35], the concentration that the cells were exposed to was determined to be ∼13.5 µM, which was the concentration used in the simulations. The concentration change in ACA (or cyanide) was simulated with a step-function for an external signal in the model using Wolfram Mathematica v.11.3 (Wolfram Research, Inc., Champaign, IL).

Phase shift analysis

Acquired images were processed and analysed as described previously [13]. The instantaneous phases of the oscillations of the individual cells, φ(t), were extracted from the NADH signals by means of a Hilbert transform of the signals in MATLAB (The MathWorks Inc.) [31,36]. The phases are defined in the interval [−π, π] and a phase of 0 represents a maximum in the concentration of NADH of the unperturbed signal. The oscillation frequency of the individual cells, ω, was calculated as a mean from time intervals 0–5, 5–10 and 20–25 min. For every cell, the phase shift relative to unperturbed oscillations was calculated using where the phase change under perturbation is given by and the expected unperturbed phase change by (with the time between and ). Phase shifts were determined 11 s after the start of the perturbation. Phase changes were mapped to the interval [−π, π] by adding or subtracting 2π where required. A similar procedure was followed for analysis of the simulation results using timecourse solutions of the oscillations with perturbations in Wolfram Mathematica v.11.3.

The degree of entrainment before and after the change in the simulated entraining signal was characterised using the order parameter calculated for the phases, , of NADH and external ACA (i.e. the order parameter is never one, but stability indicates locked phases). The degree of entrainment of intrinsic oscillation by an external signal is characterised by the order parameter [37]: 
formula
(1.1)

Results

Experimental observation and model simulation

Using our detailed model for glycolytic oscillations in isolated yeast cells (the Gustavsson model [20]), we generated a PRC for cyanide perturbations in a microfluidic chamber described experimentally in [31] (see Methods). Here, phase shifts were induced by switching between solutions with and without 5 mM cyanide in a microfluidic flow chamber. Since the yeast cells were not synchronised, their respective intrinsic phases were different at the time of perturbation. An experimental PRC can be constructed by aggregating the results of the individual cells, for each of which an instantaneous oscillation phase pre-, and post perturbation is determined, using its fluorescent NADH signal (data published in [31]; see Methods).

Our simulation result showed good agreement with experimental data (see Supplementary Figure S1) and indicated a region of phase advancement that was underrepresented in the data set of [31]. As described in [33], cyanide binds to ACA as well as many intracellular carbonyl-containing compounds [pyruvate and dihydroxyacetone phosphate (DHAP)]. The removal of cyanide therefore corresponds to, inter alia, an increase in ACA. This multi-target effect of cyanide, which complicates the mechanistic interpretation of the PRC, prompted us to investigate the effect of direct ACA perturbations.

Subsequently, direct ACA perturbations were applied to individually oscillating yeast cells, via media switching in a microfluidic chamber (see Methods). Phase shifts for the individual cells are shown as data points in Figure 1, and although heterogeneity between the individual cells leads to a considerable scatter in the data points, there is a clear dependence of the phase shift on the phase of the cells at the time of perturbation (where a phase of 0 corresponds to a maximum in NADH in the unperturbed signal). A small phase shift was observed if cells have a phase of 0 to π/2 before the perturbation, a strong positive phase shift (advancement) between 0 and −π/2, and a strong negative shift (retardation) from −π/2 to −π and from π/2 to π. Qualitatively, we see similar behaviour to the cell suspension experimental results in [30], with phase shifts most pronounced around the inflection point of rising NADH.

Single cell phase shifts in response to ACA perturbations.

Figure 1.
Single cell phase shifts in response to ACA perturbations.

The shift in the instantaneous phase of NADH is plotted as a function of the NADH phase where the ACA perturbation was applied, with a phase of 0 defined at the NADH concentration peak in the unperturbed signal. The experimental data points are compiled from single cells subjected to ACA perturbations. The grey shaded area denotes the solution space for Monte-Carlo simulations of the Gustavsson model with a 10% variation in model parameters. The red curve indicates the phase response of the model with the original parameter values (type 1) and the blue curve the response for a parameter set in the Monte-Carlo simulation leading to a type 0 response.

Figure 1.
Single cell phase shifts in response to ACA perturbations.

The shift in the instantaneous phase of NADH is plotted as a function of the NADH phase where the ACA perturbation was applied, with a phase of 0 defined at the NADH concentration peak in the unperturbed signal. The experimental data points are compiled from single cells subjected to ACA perturbations. The grey shaded area denotes the solution space for Monte-Carlo simulations of the Gustavsson model with a 10% variation in model parameters. The red curve indicates the phase response of the model with the original parameter values (type 1) and the blue curve the response for a parameter set in the Monte-Carlo simulation leading to a type 0 response.

We simulated the phase-dependent effect of the external ACA perturbation on the NADH phase shift in our mathematical model, using a similar analysis as was used for the experimental data. The model parameters are based on experimentally determined values (i.e. mean values for the population), and to simulate the heterogeneity between the individual cells, we allowed for a 10% variability in model parameters. A Monte-Carlo simulation resulted in the grey shaded area depicted in Figure 1 for 1000 parameter sets (at this point, the addition of 100 more sets resulted in only a 1% change in the area). The excellent prediction of the existing mathematical model, shown as a red curve in Figure 1 for the original parameter values, and with a grey shaded area when taking cellular heterogeneity into account, indicates that the mathematical model captures the system sufficiently well to warrant a detailed model analysis to unravel the mechanism underlying the phase response. Note that PRCs showing type 1 (as obtained with the default parameter set, shown in red) and type 0 (a typical example shown in blue) phase shift responses (as defined by Winfree [1]), fall within the grey shaded solution space of the model in Figure 1.

Transduction of the perturbation

To analyse the effect of a perturbation on the glycolytic network, a simulation for a short ACA pulse (<1 s, to minimise the effect of the network relaxation on the absorption of the perturbation) was run in the Gustavsson model. The final network phase shift was determined after the perturbation had been propagated through the network and the limit cycle had been restored. The result of the short ACA pulse was very similar to what was observed in Figure 1, showing both type-1 and type-0 PRCs, depending on the magnitude of the perturbation (PRCs are shown as green and blue curves in Supplementary Figure S2 and the perturbations that lead to the two types similarly in Supplementary Figure S3). In addition, we show in parametric plots how the phase at the time of perturbation can lead to either phase delay or phase advancement (Supplementary Figure S3).

The propagation of the ACA perturbation into the network was analysed by following the subsequent changes in internal metabolite concentrations and reaction rates. As depicted schematically in Figure 2 (see Supplementary Figure S5 for time courses), the rise in external ACA causes a rapid rise in internal ACA, which accelerates the alcohol dehydrogenase reaction (ADH) causing an increased oxidation rate of NADH, and a decrease in its concentration. This decrease in NADH (and concomitant increase in NAD) concentration stimulates the glyceraldehyde 3-phosphate dehydrogenase (GAPDH) reaction (as also proposed in [30]). The activity of GAPDH is closely followed by the phosphoglycerate kinase (PGK) reaction, leading to an increased ATP concentration. The translation of the ACA perturbation to internal ATP concentration is a rapid signal transduction step with no significant delay, running mostly via the cofactors [38,39].

Signal transduction of an external ACA perturbation in the network.

Figure 2.
Signal transduction of an external ACA perturbation in the network.

An increased external ACA enters rapidly into the cell (step 1) and activates ADH, thereby decreasing NADH and increasing NAD. This leads to an acceleration of GAPDH and the tightly coupled PGK, leading to an increase in ATP. The signal transduction of ACA to the cofactors NADH and ATP (step 2) is fast. Step 3, transduction of the ATP signal to the kinases and subsequent relaxation in the pathway is slow, and ultimately leads to the phase shift response via its effect on PFK.

Figure 2.
Signal transduction of an external ACA perturbation in the network.

An increased external ACA enters rapidly into the cell (step 1) and activates ADH, thereby decreasing NADH and increasing NAD. This leads to an acceleration of GAPDH and the tightly coupled PGK, leading to an increase in ATP. The signal transduction of ACA to the cofactors NADH and ATP (step 2) is fast. Step 3, transduction of the ATP signal to the kinases and subsequent relaxation in the pathway is slow, and ultimately leads to the phase shift response via its effect on PFK.

The increased ATP concentration affects the other kinases in the system. The glucokinase (GLK) activity is accelerated (ATP is a substrate), and although pyruvate kinase (PYK) activity is modestly inhibited by ATP (a product), this inhibition is overcome by an increase in PEP (a substrate) resulting from the initial increase in GAPDH/PGK activities. The adenylate kinase (AK) operates close to equilibrium, and the increased ATP concentration leads to a decrease in AMP concentration. PFK is sensitive for changes in ATP which acts both as a substrate and as an allosteric inhibitor, and also to changes in AMP which acts as an activator of the enzyme. The rise in ATP concentration and the decrease in the AMP concentration has the net effect of inhibiting the PFK reaction (see Supplementary Figure S6). We subsequently compared the immediate local changes in kinase activities with the final system phase changes and noted that the local changes in GLK and PFK activities were close to the phase shift in the relaxed system, while there was no such correlation for the PYK (see Supplementary Figures S6, S7). To finally distinguish between the contributions of the GLK and PFK we made a perturbation analysis, while keeping the respective enzymes at their original limit cycle activities for up to half a period and then releasing them again. Doing so for the GLK had no effect on the system's phase shift, while clamping the PFK had a significant effect. From this, we could conclude that the PFK drives the phase shift, and that GLK simply follows the phase shift dictated by the PFK, i.e. after releasing the GLK from its clamped state it immediately jumps to the system's phase. The method we used for the mechanistic interpretation of the phase shift response is generally applicable if a detailed mathematical model is available.

Phase shift response to synchronisation

PRCs are important tools for studying synchronisation between oscillators. For entrainment or synchronisation, a phase shift should show a qualitatively correct adaptation to the perturbation, i.e. a phase advancement when the perturbation phase is advanced. This criterion can be readily applied to excitation oscillators such as nervous cells, where a perturbation spike can be related to the excitation of the system. For a more sinusoidal oscillation, it is not possible to define a phase for a single perturbation spike, and we therefore compared the pulse-induced phase shifts (as determined in Supplementary Figure S2 for positive pulses and Supplementary Figure S8 for negative pulses) with continuous external signals that were either delayed or advanced in phase relative to the intrinsic ACA. In Figure 3, we plot the NADH and ACA concentrations for the non-perturbed oscillation (blue and orange curves), and the ACA concentration for an entraining signal with an advanced phase (green curve) or a retarded phase (red curve). In addition, we show in grey the areas where the phase shift is positive (phase advancement) upon a positive (Figure 3A,B) or negative (Figure 3C,D) ACA pulse. These plots show that for the largest part of the period the cells adapt correctly, i.e. they show a phase advancement when ACA is perturbed in the direction of the advanced signal. For example, when the advanced signal has a higher ACA concentration (green curve above orange curve), the cells show phase advancement as they would for positive ACA pulses, indicated by the upward green arrows in the grey area in Figure 3A. When the advanced oscillator has a lower ACA concentration, the cells show phase advancement as they would for negative ACA pulses indicated by the downward green arrows in the grey area in Figure 3C. The same holds for the delayed signal shown in red in Figure 3B,D.

The direction of the phase shift response relative to the entraining signal.

Figure 3.
The direction of the phase shift response relative to the entraining signal.

The blue and orange curves in (A)–(D) indicate the NADH and ACA concentrations for the non-perturbed Gustavsson model. (A and C) The effect of an advanced external ACA signal (green curve) is displayed. The green arrows show that the advanced signal reflects an ACA pulse (arrows pointing upwards), or an ACA removal (negative pulse, indicated by arrows pointing downwards). The grey regions indicate the areas where an ACA pulse (A and B), or an ACA removal (C and D), leads to phase advancement, as determined with a single pulse perturbation, up (A and B) or down (C and D). The white regions indicate the areas where an ACA pulse (A and B), or an ACA removal (C and D), leads to phase delay. The left column shows that the regions of increased and decreased ACA, coming from an advanced signal, predominantly lead to phase advancement (mostly in grey areas) and the cell's phase is advanced to adapt to the new phase. The right column shows that the regions of increased and decreased ACA, coming from a delayed signal, predominantly lead to phase delay (white regions) and the cell's phase is retarded to adapt to the new phase of the entraining signal.

Figure 3.
The direction of the phase shift response relative to the entraining signal.

The blue and orange curves in (A)–(D) indicate the NADH and ACA concentrations for the non-perturbed Gustavsson model. (A and C) The effect of an advanced external ACA signal (green curve) is displayed. The green arrows show that the advanced signal reflects an ACA pulse (arrows pointing upwards), or an ACA removal (negative pulse, indicated by arrows pointing downwards). The grey regions indicate the areas where an ACA pulse (A and B), or an ACA removal (C and D), leads to phase advancement, as determined with a single pulse perturbation, up (A and B) or down (C and D). The white regions indicate the areas where an ACA pulse (A and B), or an ACA removal (C and D), leads to phase delay. The left column shows that the regions of increased and decreased ACA, coming from an advanced signal, predominantly lead to phase advancement (mostly in grey areas) and the cell's phase is advanced to adapt to the new phase. The right column shows that the regions of increased and decreased ACA, coming from a delayed signal, predominantly lead to phase delay (white regions) and the cell's phase is retarded to adapt to the new phase of the entraining signal.

Control of adaptation time

In addition to phase responses, we also considered the adaptation time to reach a new limit cycle after the system is exposed to an external entrainment signal. The time courses for the ACA perturbations show a fast initial effect on the cofactor levels and a subsequent slow relaxation back to the limit cycle. We analysed the importance of the kinases on the adaptation time to an external entraining signal. We simulated an external ACA entrainment signal (2 µM, sinusoidal with the same period as the model), and after phase locking, changed the signal phase by 180°. As can be seen in Figure 4, it takes considerable time for the wild-type cells to adapt themselves to the new signal after the phase flip at time 0. This wild-type response is indicated in terms of amplitude and phase restoration separately in Figure 4A (purple and blue curves, respectively). Here, the order parameter was used as a measure of phase locking between NADH and ACA (see Methods), i.e. a constant order parameter indicates a stable lock. The ratio is used as a measure of amplitude entrainment. As also shown in [15], the glycolytic oscillator shows more rapid amplitude than phase restoration when the entraining signal's phase is flipped by 180°. The same can be seen when considering the blue trajectory on the amplitude and phase-polar plot in Figure 4B. The time interval between successive dots is one period of the oscillation and the outer circle indicates fully restored amplitude. A fully restored phase entrainment is attained at 271° (with the ACA signal leading). The ability of PFK, GLK and PYK to change the adaptation time of the cell to the new entraining signal was tested by, in turn, briefly (0.05 × period ≈ 2.5 s) changing their specific activity by 10% at the point where the entraining signal is changed, and following the network's subsequent return to the limit cycle. This analysis was repeated across one period of oscillations and overall a change in the activity of PFK induced the largest change in the relaxation time to the limit cycle. This was also verified for a perturbation time of 0.01 × period (∼0.5 s). The effect of the concomitant PFK perturbation on amplitude and phase is depicted as the red and orange curves in Figure 4A and the orange trajectory in Figure 4B. Both amplitude and phase show faster adaptation to entrainment. Figure 4C shows the adaptation to the entrained limit cycle as reflected in the timecourse of PFK and its substrate and product (F6P and F16bP). The green curve depicts the timecourse in the case where there is no 180° phase flip in the entraining signal. The blue curve depicts the wild-type adaptation to the new signal and the orange curve the adaptation in the case where PFK is briefly perturbed. Clearly, the orange and green curves have settled at a 180° phase difference long before the wild-type adaptation (blue) manages to adapt to the new entraining signal and overlaps with the orange curve.

PFK can shorten the adaptation time to a change in an external entraining signal.

Figure 4.
PFK can shorten the adaptation time to a change in an external entraining signal.

At time = 0, a 180° phase shift in an external entraining ACA signal is made and the effect of a PFK perturbation on the adaptation followed in the Gustavsson model. (A) Smoothed phase order parameter without (blue) and with (orange) an additional PFK perturbation. The amplitude recovery is shown without (purple) and with (red) the same perturbation. The PFK perturbation shortens the adaptation time of the amplitude and the phase. (B) Phase and amplitude response on a polar plot where the angle indicates the phase difference between the forcing ACA signal and the NADH signal (with the ACA signal leading), and the distance from the centre indicates the NADH amplitude relative to the fully entrained NADH amplitude with the outer ring being 1. Here, arrows indicate the direction of movement and the time between coloured dots is 1 period. Here, it is also evident that the PFK perturbation (orange trajectory) leads to a faster adaptation to the entraining signal. (C) Effect of the phase shift on the timecourses of PFK activity, and F6P and F16bP concentrations. The green curve always indicates the timecourse without a phase shift, the blue curve the timecourse with a phase shift and a subsequent natural adaptation (right panel), and the orange curve the timecourse with a phase shift and a PFK perturbation leading to a faster adaptation (left panel) to be 180° changed relative to the green curve. Arrows indicate the direction of phase advancement induced by the PFK perturbation (orange) and the natural adaptation (blue) which manifests as a slow phase retardation.

Figure 4.
PFK can shorten the adaptation time to a change in an external entraining signal.

At time = 0, a 180° phase shift in an external entraining ACA signal is made and the effect of a PFK perturbation on the adaptation followed in the Gustavsson model. (A) Smoothed phase order parameter without (blue) and with (orange) an additional PFK perturbation. The amplitude recovery is shown without (purple) and with (red) the same perturbation. The PFK perturbation shortens the adaptation time of the amplitude and the phase. (B) Phase and amplitude response on a polar plot where the angle indicates the phase difference between the forcing ACA signal and the NADH signal (with the ACA signal leading), and the distance from the centre indicates the NADH amplitude relative to the fully entrained NADH amplitude with the outer ring being 1. Here, arrows indicate the direction of movement and the time between coloured dots is 1 period. Here, it is also evident that the PFK perturbation (orange trajectory) leads to a faster adaptation to the entraining signal. (C) Effect of the phase shift on the timecourses of PFK activity, and F6P and F16bP concentrations. The green curve always indicates the timecourse without a phase shift, the blue curve the timecourse with a phase shift and a subsequent natural adaptation (right panel), and the orange curve the timecourse with a phase shift and a PFK perturbation leading to a faster adaptation (left panel) to be 180° changed relative to the green curve. Arrows indicate the direction of phase advancement induced by the PFK perturbation (orange) and the natural adaptation (blue) which manifests as a slow phase retardation.

Discussion

Phase-response curves provide a powerful and experimentally assessable tool to study autonomous oscillators. They give important insight into overall systems dynamics, but the advantage that the analysis can be performed on the intact system often limits the number of perturbations that can be made, and the interpretation of the curves often lacks mechanistic understanding. The strongest applications of phase-response curves are arguably in the field of neuroscience where strong mathematical models underlie the perturbation experiments. Many mathematical models exist for yeast glycolytic oscillations, although most of them are core models, and do not capture the full biochemical mechanisms underlying the oscillations. PRCs have not been used extensively to study yeast glycolytic oscillations in intact cells (the analyses of which were limited to synchronised populations studies for a long time), but have been studied in cell-free extracts, and analysed with mathematical models. The cell-free extracts studies have a clear advantage that more direct perturbations can be made, i.e. not only external variables as in the case of intact cells.

In [2123], it was shown that in oscillating yeast cell extracts the addition of any of the adenosine phosphates causes a phase shift varying in direction and magnitude depending on where in the phase the addition is made. Specifically, it was seen that if ADP or AMP is added at a minimum of NADH, a maximal phase shift is observed, the equivalence being a result of the interconversion by AK. For ATP alone this effect was reversed. The interconversion of ADP to AMP leads to ATP formation as well, but the AMP effect seems to overrule that of ATP. In [23], the adenosine phosphate effects were postulated to be due to AMP being an allosteric activator, and ATP an allosteric inhibitor of PFK, in addition to their roles as substrate activators (ADP) and product inhibitors (ATP) of the kinases in the lower part of glycolysis.

In [28], Golbeter found his model of PFK allosterically activated by its product, ADP, to be sufficient for describing the observed phase shifts when adding ADP. Addition of the stimulating product when the enzyme is at maximal activity results in minimal phase shift since no significant additional activation is possible. When ADP is added near the minimal activity of PFK, maximal phase shift is observed. In this case, the premature acceleration of the reaction leads to a premature decline in the substrate (ATP) which has to re-accumulate before synthesis of the product can again occur leading to phase retardation (see Figure 5A).

Comparison of the PRCs in the Goldbeter and Gustavsson models.

Figure 5.
Comparison of the PRCs in the Goldbeter and Gustavsson models.

A pulse of the allosteric activator of PFK (ADP in the case of the Goldbeter model [28] (A) and AMP in the case of the Gustavsson model (B)) is administered to induce a phase response (red curves). Timecourses for PFK and ADP (or AMP) are shown for one period as an overlay. The zero phase at perturbation is here determined by the peak in PFK activity.

Figure 5.
Comparison of the PRCs in the Goldbeter and Gustavsson models.

A pulse of the allosteric activator of PFK (ADP in the case of the Goldbeter model [28] (A) and AMP in the case of the Gustavsson model (B)) is administered to induce a phase response (red curves). Timecourses for PFK and ADP (or AMP) are shown for one period as an overlay. The zero phase at perturbation is here determined by the peak in PFK activity.

In the Gustavsson single cell model analysed here, PFK is activated by its substrate, fructose-6-phosphate, in addition to being allosterically activated by AMP, formed via the action of AK from the PFK product ADP. In addition, PFK is inhibited by its co-substrate, ATP, in the concentration range that spans the ATP oscillation amplitude. This leads to a more complex regulation pattern; however, we did test the phase response in the Gustavsson model upon an AMP pulse (Figure 5B). On comparing with Figure 5A, it is clear that the Gustavsson model shows a similar response as a function of PFK activity. At minimum activity, a pulse of the allosteric activator leads to phase retardation and at maximum activity, the phase shift is much reduced.

In contrast with the PFK acceleration by an ADP pulse as described in [28] (or an AMP pulse as examined above), the effect of the ACA pulse translates into the removal of AMP with a concomitant increase in ATP. The net effect is an instantaneous inhibition of PFK despite an acceleration of the upper and lower parts of glycolysis. We see minimal phase shift if PFK is inhibited near its minimum activity and maximal phase shift if this happens around maximal activity. If the activity of PFK is increasing at the time (to the left of the peak), this inhibition leads to phase retardation. Around maximal activity, this inhibition leads to premature deceleration and phase advancement (see Supplementary Figure S9). In the region of decreasing activity to the far right of the peak, the activity is inhibited to such an extent that it overshoots the minimum leading also to phase retardation. We also tested our model for the similarity between an AMP addition and a short ACA sink (i.e. negative pulse) and obtained good agreement (see Supplementary Figure S10).

As noted in [7], PRCs for many biological systems reflect a degree of similarity in phase advancement and retardation regions, and entrainment properties, which implies that these curves, per se, give no information about the detailed mechanisms underlying the phase shifts and synchronisation. As we have demonstrated, analysis of the PRC using a detailed mechanistic model of the system can give information on features of the system required for the observed phase shifts and by implication also required for individual oscillators to exhibit synchronisation behaviour. The kinetic properties of PFK, which have been shown to be necessary for oscillations to occur [20], also play a key role in the phase response of the glycolytic network, and the adaptation time to an entraining signal.

One of the autonomous oscillatory systems for which effects of perturbations are widely known, is the circadian clock, with shift work and travelling through time zones as frequently occurring perturbations with annoying side effects. PRCs of the circadian clock have been investigated experimentally by demonstrating the effect that light pulses have on the human circadian phase [4045], and have been studied with mathematical models [10]. These theoretical models have a molecular basis and incorporate a negative autoregulation of gene expression [46] and can describe the experimental data. Whereas these models have been suggested to be useful to analyse disorders of circadian rhythms, as far as we know they have not been used for fine-tuning of small phase shifts, resulting from travelling over time zones. In particular, it would be important that in addition to the molecular components in the model, the perturbation to be made, e.g. melatonin intake, should be part of the model. The initiative toward a detailed computational model for the mammalian circadian clock [10] will not only provide a fundamental understanding of the important components but also pinpoint targets for possible phase intervention.

Abbreviations

     
  • ACA

    acetaldehyde

  •  
  • ADH

    alcohol dehydrogenase

  •  
  • AK

    adenylate kinase

  •  
  • DHAP

    dihydroxyacetone phosphate

  •  
  • GAPDH

    glyceraldehyde 3-phosphate dehydrogenase

  •  
  • GLK

    glucokinase

  •  
  • PFK

    Phosphofructokinase

  •  
  • PGK

    phosphoglycerate kinase

  •  
  • PRC

    phase response curve

  •  
  • PYK

    pyruvate kinase

Author Contribution

A.-K.G. and M.M.-B. performed the experiments and analysed the data and under the guidance of C.B.A. and M.G., D.D.v.N. and J.L.S. performed model simulations, analyses and interpretations. D.D.v.N. and J.L.S. wrote the manuscript. All authors discussed and commented on the manuscript.

Funding

We acknowledge financial assistance from the DST/NRF in South Africa, particularly for funding the SARChI initiative (NRF-SARCHI-82813) and for Thuthuka funding (TTK14051967526). This work was furthermore funded by the Swedish Research Council and the Carl Tryggers Foundation for Scientific Research.

Competing Interests

The Authors declare that there are no competing interests associated with the manuscript.

References

References
1
Winfree
,
A.T.
(
1970
)
Integrated view of resetting a circadian clock
.
J. Theor. Biol.
28
,
327
374
2
Winfree
,
A.T.
(
2001
) The geometry of biological time. In
Interdisciplinary Applied Mathematics
,
2nd edn
, vol.
12
,
Springer-Verlag
,
New York, NY
3
Pittendrigh
,
C.S.
and
Daan
,
S.
(
1976
)
A functional analysis of circadian pacemakers in nocturnal rodents - I. The stability and lability of spontaneous frequency
.
J. Comp. Physiol. A
106
,
223
252
4
Reyes
,
A.D.
and
Fetz
,
E.E.
(
1993
)
Effects of transient depolarizing potentials on the firing rate of cat neocortical neurons
.
J. Neurophysiol.
69
,
1673
1683
5
Guevara
,
M.R.
,
Shrier
,
A.
and
Glass
,
L.
(
1986
)
Phase resetting of spontaneously beating embryonic ventricular heart cell aggregates
.
Am. J. Physiol.
251
,
H1298
H1305
6
Schultheiss
,
N.W.
,
Prinz
,
A.A.
and
Butera
,
R.J.
(
2011
)
Phase Response Curves in Neuroscience. Theory, Experiment, and Analysis
,
Springer
,
New York
7
Granada
,
A.
,
Hennig
,
R.M.
,
Ronacher
,
B.
,
Kramer
,
A.
and
Herzel
,
H.
(
2009
)
Phase response curves elucidating the dynamics of coupled oscillators
.
Methods Enzymol.
454
,
1
27
8
Saunders
,
D.S.
and
Thomson
,
E.J.
(
1977
)
‘Strong’ phase response curve for the circadian rhythm of locomotor activity in a cockroach (Nauphoeta cinerea)
.
Nature
270
,
241
243
9
Burgess
,
H.J.
,
Revell
,
V.L.
,
Molina
,
T.A.
and
Eastman
,
C.I.
(
2010
)
Human phase response curves to three days of daily melatonin: 0.5 mg Versus 3.0 mg
.
J. Clin. Endocrinol. Metab.
95
,
3325
3331
10
Leloup
,
J.C.
and
Goldbeter
,
A.
(
2003
)
Toward a detailed computational model for the mammalian circadian clock
.
Proc. Natl Acad. Sci. U.S.A.
100
,
7051
7056
11
Kralemann
,
B.
,
Frühwirth
,
M.
,
Pikovsky
,
A.
,
Rosenblum
,
M.
,
Kenner
,
T.
,
Schaefer
,
J.
et al.  (
2013
)
In vivo cardiac phase response curve elucidates human respiratory heart rate variability
.
Nat. Commun.
4
,
2418
12
Richard
,
P.
(
2003
)
The rhythm of yeast
.
FEMS Microbiol. Rev.
27
,
547
557
13
Gustavsson
,
A.K.
,
van Niekerk
,
D.D.
,
Adiels
,
C.B.
,
du Preez
,
F.B.
,
Goksör
,
M.
and
Snoep
,
J.L.
(
2012
)
Sustained glycolytic oscillations in individual isolated yeast cells
.
FEBS J.
279
,
2837
2847
14
Richard
,
P.
,
Bakker
,
B.M.
,
Teusink
,
B.
,
van Dam
,
K.
and
Westerhoff
,
H.
(
1996
)
Acetaldehyde mediates the synchronization of sustained glycolytic oscillations in populations of yeast cells
.
Eur. J. Biochem.
235
,
238
241
15
Danø
,
S.
,
Madsen
,
M.F.
and
Sørensen
,
P.G.
(
2007
)
Quantitative characterization of cell synchronization in yeast
.
Proc. Natl Acad. Sci. U.S.A.
104
,
12732
12736
16
Sel'kov
,
E.E.
(
1975
)
Stabilization of energy charge, generation of oscillations and multiple steady states in energy metabolism as a result of purely stoichiometric regulation
.
Eur. J. Biochem.
59
,
151
157
17
Hess
,
B.
(
1979
)
The glycolytic oscillator
.
J. Exp. Biol.
81
,
7
14
PMID:
[PubMed]
18
du Preez
,
F.B.
,
van Niekerk
,
D.D.
,
Kooi
,
B.
,
Rohwer
,
J.M.
and
Snoep
,
J.L.
(
2012
)
From steady-state to synchronized yeast glycolytic oscillations I: model construction
.
FEBS J.
279
,
2810
2822
19
du Preez
,
F.B.
,
van Niekerk
,
D.D.
and
Snoep
,
J.L.
(
2012
)
From steady-state to synchronized yeast glycolytic oscillations II: model validation
.
FEBS J.
279
,
2823
2836
20
Gustavsson
,
A.K.
,
van Niekerk
,
D.D.
,
Adiels
,
C.B.
,
Kooi
,
B.
,
Goksör
,
M.
and
Snoep
,
J.L.
(
2014
)
Allosteric regulation of phosphofructokinase controls the emergence of glycolytic oscillations in isolated yeast cells
.
FEBS J.
281
,
2784
2793
21
Chance
,
B.
,
Schoener
,
B.
and
Elsaesser
,
S.
(
1964
)
Control of the waveform of oscillations of the reduced pyridine nucleotide level in a cell-free extract
.
Proc. Natl Acad. Sci. U.S.A.
52
,
337
341
22
Chance
,
B.
,
Schoener
,
B.
and
Elsaesser
,
S.
(
1965
)
Metabolic control phenomena involved in damped sinusoidal oscillations of reduced diphosphopyridine nucleotide in a cell-free extract of Saccharomyces carlsbergensis
.
J. Biol. Chem.
240
,
3170
3181
PMID:
[PubMed]
23
Hess
,
B.
,
Boiteux
,
A.
and
Krüger
,
J.
(
1969
)
Cooperation of glycolytic enzymes
.
Adv. Enzyme Regul.
7
,
149
167
24
Pye
,
E.K.
(
1973
) Glycolytic oscillations in cells and extracts of yeast-some unsolved problems. In
Biological and Biochemical Oscillators
, pp.
269
284
,
Academic Press
25
Goldbeter
,
A.
and
Lefever
,
R.
(
1972
)
Dissipative structures for an allosteric model: application to glycolytic oscillations
.
Biophys. J.
12
,
1302
1315
26
Goldbeter
,
A.
(
1974
)
Modulation of the adenylate energy charge by sustained metabolic oscillations
.
FEBS Lett.
43
,
327
330
27
Boiteux
,
A.
,
Goldbeter
,
A.
and
Hess
,
B.
(
1975
)
Control of oscillating glycolysis of yeast by stochastic, periodic, and steady source of substrate: a model and experimental study
.
Proc. Natl Acad. Sci. U.S.A.
72
,
3829
3833
28
Goldbeter
,
A.
(
1997
)
Biochemical Oscillations and Cellular Rhythms
,
Cambridge University Press
,
Cambridge, U.K.
29
Madsen
,
M.F.
,
Danø
,
S.
and
Sørensen
,
P.G.
(
2005
)
On the mechanisms of glycolytic oscillations in yeast
.
FEBS J.
272
,
2648
2660
30
Betz
,
A.
and
Becker
,
J.U.
(
1975
)
Phase dependent phase shifts induced by pyruvate and acetaldehyde in oscillating NADH of yeast cells
.
J. Interdiscip. Cycle Res.
6
,
167
173
31
Gustavsson
,
A.K.
,
Adiels
,
C.B.
,
Mehlig
,
B.
and
Goksör
,
M.
(
2015
)
Entrainment of heterogeneous glycolytic oscillations in single cells
.
Sci. Rep.
5
,
9404
9407
32
Yates
,
W.F.
and
Heider
,
R.L.
(
1952
)
The dissociation of lactonitrile in aqueous solution
.
J. Am. Chem. Soc.
74
,
4153
4155
33
Hald
,
B.O.
,
Nielsen
,
A.G.
,
Tortzen
,
C.
and
Sørensen
,
P.G.
(
2015
)
Cyanohydrin reactions enhance glycolytic oscillations in yeast
.
Biophys. Chem.
200–201
,
18
26
34
Richard
,
P.
,
Teusink
,
B.
,
Westerhoff
,
H.
and
van Dam
,
K.
(
1993
)
Around the growth phase transition S. cerevisiae's make-up favours sustained oscillations of intracellular metabolites
.
FEBS Lett.
318
,
80
82
35
Hald
,
B.O.
,
Smrcinova
,
M.
and
Sørensen
,
P.G.
(
2012
)
Influence of cyanide on diauxic oscillations in yeast
.
FEBS J.
279
,
4410
4420
36
Rosenblum
,
M.
,
Pikovsky
,
A.
and
Kurths
,
J.
(
1996
)
Phase synchronization of chaotic oscillators
.
Phys. Rev. Lett.
76
,
1804
1807
37
Shinomoto
,
S.
and
Kuramoto
,
Y.
(
1986
)
Phase transitions in active rotator systems
.
Prog. Theor. Phys.
75
,
1105
1110
38
Wolf
,
J.
,
Passarge
,
J.
,
Somsen
,
O.J.
,
Snoep
,
J.L.
,
Heinrich
,
R.
and
Westerhoff
,
H.
(
2000
)
Transduction of intracellular and intercellular dynamics in yeast glycolytic oscillations
.
Biophys. J.
78
,
1145
1153
39
Bier
,
M.
,
Bakker
,
B.M.
and
Westerhoff
,
H.
(
2000
)
How yeast cells synchronize their glycolytic oscillations: a perturbation analytic treatment
.
Biophys. J.
78
,
1087
1093
40
Honma
,
K.
and
Honma
,
S.
(
1988
)
A human phase response curve for bright light pulses
.
Jpn J. Psychiatry Neurol.
42
,
167
168
41
Czeisler
,
C.A.
,
Kronauer
,
R.E.
,
Allan
,
J.S.
,
Duffy
,
J.F.
,
Jewett
,
M.E.
,
Brown
,
E.N.
et al.  (
1989
)
Bright light induction of strong (type 0) resetting of the human circadian pacemaker
.
Science
244
,
1328
1333
42
Minors
,
D.S.
,
Waterhouse
,
J.M.
and
Wirz-Justice
,
A.
(
1991
)
A human phase-response curve to light
.
Neurosci. Lett.
133
,
36
40
43
Zeitzer
,
J.M.
,
Dijk
,
D.J.
,
Kronauer
,
R.E.
,
Brown
,
E.N.
and
Czeisler
,
C.A.
(
2000
)
Sensitivity of the human circadian pacemaker to nocturnal light: melatonin phase resetting and suppression
.
J. Physiol.
526
,
695
702
44
Khalsa
,
S.B.S.
,
Jewett
,
M.E.
,
Cajochen
,
C.
and
Czeisler
,
C.A.
(
2004
)
A phase response curve to single bright light pulses in human subjects
.
J. Physiol.
549
,
945
952
45
St Hilaire
,
M.A.
,
Gooley
,
J.J.
,
Khalsa
,
S.B.S.
,
Kronauer
,
R.E.
,
Czeisler
,
C.A.
and
Lockley
,
S.W.
(
2012
)
Human phase response curve to a 1 h pulse of bright white light
.
J. Physiol.
590
,
3035
3045
46
Leloup
,
J.C.
and
Goldbeter
,
A.
(
1998
)
A model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins
.
J. Biol. Rhythms
13
,
70
87

Supplementary data