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)
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 [1–6], but it can also be used to quantify coupling strengths and to predict complex dynamics resulting from perturbations . 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 , 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 , 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 , 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  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 ; oscillations in single cells ; 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  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,18–20]. 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 [21–23] and ACA [22,24]. In [25–28], 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 .
Since it has been shown that single cells in isolation do oscillate , 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 , 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  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 , 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 , 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 , 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  to analyse the cyanide perturbations published in . Although cyanide does trap ACA via lactonitrile formation , it binds also to many other metabolites . 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.
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 . 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.
The experimental setup and procedure were described in detail in . 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 . 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  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 , 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 . 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.
Experimental observation and model simulation
Using our detailed model for glycolytic oscillations in isolated yeast cells (the Gustavsson model ), we generated a PRC for cyanide perturbations in a microfluidic chamber described experimentally in  (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 ; 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 . As described in , 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 , with phase shifts most pronounced around the inflection point of rising NADH.
Single cell phase shifts in response to ACA perturbations.
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 ), 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 ). 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.
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.
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 , 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.
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 [21–23], 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 , 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 , 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.
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  (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 , 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 , 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 [40–45], and have been studied with mathematical models . These theoretical models have a molecular basis and incorporate a negative autoregulation of gene expression  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  will not only provide a fundamental understanding of the important components but also pinpoint targets for possible phase intervention.
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.
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.
The Authors declare that there are no competing interests associated with the manuscript.