Emergence of collective oscillations in adaptive cells

Collective oscillations of cells in a population appear under diverse biological contexts. Here, we establish a set of common principles by categorising the response of individual cells against a time-varying signal. A positive intracellular signal relay of sufficient gain from participating cells is required to sustain the oscillations, together with phase matching. The two conditions yield quantitative predictions for the onset cell density and frequency in terms of measured single-cell and signal response functions. Through mathematical constructions, we show that cells that adapt to a constant stimulus fulfil the phase requirement by developing a leading phase in an active frequency window that enables cell-to-signal energy flow. Analysis of dynamical quorum sensing in several cellular systems with increasing biological complexity reaffirms the pivotal role of adaptation in powering oscillations in an otherwise dissipative cell-to-cell communication channel. The physical conditions identified also apply to synthetic oscillatory systems.

. Letters in blue denote metabolites, while those in red are the reactions. Directional (bidirectional) arrows indicate irreversible (reversible) reactions. Abbreviations: Glco, glucose; ACE, acetaldehyde, ADH, alcohol dehydrogenase; AK, adenylate kinase; ALD, fructose-1,6-bisphosphate aldolase; BPG, 1,3-bis-phosphoglycerate; ENO, phosphopyruvate hydratase; F16P, fructose-1,6bisphosphate; F6P, fructose 6-phosphate; GAPDH, D-glyceraldehyde-3-phosphate dehydrogenase (phosphorylating); G3P, glycerol 3-phosphate; G3PDH, glycerol 3-phos-phate dehydrogenase; G6P, glucose 6-phosphate; GLYCO, glycogen branch; GLK, glucokinase (a hexokinase); P2G, 2-phosphoglycerate; P3G, 3-phosphoglycerate; PEP, phosphoenolpyruvate; PDC, pyruvate decarboxylase; PGI, glucose-6-phosphate isomerase; PFK, 6-phosphofructokinase; PGK, phosphoglycerate kinase; PGM, phosphoglycerate mutase; PYK, pyruvate kinase; PYR, pyruvate; Treha, trehalose branch; SUC, succinate branch; GLYO, glyoxylate shunt. Consider the temporal variation a t of an intracellular observable a induced by a sinusoidal signal s t from the environment at frequency ω. The amplitude of s t is assumed to be small, so that it can be treated as a perturbation. The observable a is either directly or indirectly coupled to the signal s. On the scale of a single cell, both a t and s t may contain stochastic components. In the following, we shall examine the noise averaged response of a to the deterministic part of s, i.e., the signal. As usual, we use · to denote the noise average. The following convention on forward and inverse Fourier transforms is adopted, In a steady-state, the ratio of the Fourier amplitudes ã(ω) and s(ω) defines the response functionR a (ω) ≡ ã(ω) / s(ω) , which can be separated into its realR a and imaginaryR a parts. (For a cellular variable a that follows stochastic dynamics, the ensemble averaged response is considered.) The well-known Kramers-Krönig relation from causality requirement on the response function states 2 : In the case of a perfectly adapting a, the response vanishes under a sufficiently slow stimulus, i.e., lim ω→0Ra (ω) = 0. Equation (2) then requires, Consequently,R a (ω) must change sign at least once along the frequency axis. Let φ a = − arg(R a ) be the phase of −R a (ω), with the minus sign introduced by convention. Positive and negative values ofR a thus translate to phase-lag (−π < φ a < 0) and phase-lead (0 < φ a < π) of a t over s t , respectively. By virtue of continuity, the sign change of R a (ω) is also expected in the partially adaptive case, provided the adaptation error R a (0) is sufficiently small.

Energy outflow from an adaptive channel
Auto-induced collective oscillations in a dissipative medium require an energy source. Below, we show that an active cell is able to output energy to a fluctuating s in the presence of an adaptive channel. The power of the output depends on the strength of the coupling as well as the amplitude and frequency of the fluctuating signal.
Consider a slightly more general form of Equation (1) in the Main Text where the contribution from cell j to the total thermodynamic force on s is given by O(a j ), which in general is nonlinear in a j . The work done on s by the cell in a time interval (0, L) can then be written as where O t ≡ O(a j (t)) and s t are both fluctuating quantities in general. We now consider a sinusoidal signal s t = s 0 + ∆s cos(ωt) with a small amplitude ∆s. To the first order in ∆s, we have Here O (0) t denotes the stochastic trajectory of O in the absence of the sinusoidal signal. As usual, the linear response function R O in the steady-state (ss) to a weak time-varying signal s t is introduced through where · denotes average over noise. The phase angle φ O (ω) ≡ − argR O (ω). Substituting expression (5) into Eq. (4) and taking the limit L → ∞, we obtain the time-averaged output power from the cell through this channel (omitting the subscript j),Ẇ Here the overline bar indicates averaging over time, and o (∆s) 2 denotes terms higher than second order in ∆s. Given the relation O t = O(a t ), adaptation of the cellular variable a to a slow-varying s also implies the adaption of O to the signal. The causality condition (2) applied to O t then requiresR O (ω) ≡ −|R O (ω)| sin φ O (ω) < 0 in a certain frequency range. Consequently, Eq. (7) predicts energy outflow from the channel under a periodic stimulation at these frequencies.
The discussion leading to Eq. (7) in the previous section can be easily extended to the energy outflow under an arbitrary signal variation s t with a power spectrumC s (ω), where ∆s sets the overall amplitude of signal variation. If the cell were in thermal equilibrium, a would respond passively to a time-varying signal with a phase-lag and dissipate the energy inflow generated by the stimulation. An adaptive cell, on the other hand, is able to output energy in the form of work when stimulated in the right frequency range. This form of energy outflow is different from the heat dissipation arising from keeping the system out of equilibrium as studied in Refs. [3][4][5] .

The Fluctuation-Dissipation Theorem
The fluctuation-dissipation theorem (FDT) is generally presented as an identity between the response function of a chosen variable to an external perturbation and the correlation function of the variable in question with the one that is conjugate to the perturbation 6 . For Markov systems which are of interest here, FDT holds when the detailed balance condition on the state-space transition rates is fulfilled. We refer the reader to Refs. 7,8 for a detailed discussion, including more rigorous definitions of various quantities of interest.
Assuming that the signal s affects the cell through coupling to a conjugate variable O which is proportional to the variable a of interest, i.e., O = c 0 a with c 0 a proportionality constant. In this case, FDT states that Here,C O (ω) = c 2 0 |ã(ω)| 2 is the power spectrum of O t which is always positive. Equation (9) contradicts (3), re-affirming that receptor adaptation cannot be realised without the presence of active processes inside the cell.
In Ref. 9 , adaptation through a 3-node incoherent feed-forward motif was considered. It was later shown that the topology even supports adaptation in an equilibrium setting 10 . The main difference between these models and the adaptive receptor model in the Main Text ( Fig. 3 and Methods) is that, in the former, s not only couples to a directly, but also to other intracellular variables. The conjugate variable O is then a combination of a and other intracellular variables. We leave a detailed investigation of this issue to future work.

Supplementary Note 2. A self-consistent scheme for frequency selection and oscillation amplitude determination
The thermodynamic analysis in the preceding section suggests the possibility of a positive feedback loop formed by a periodic signal and adaptive cells under generic conditions. Collective oscillations emerge when signal amplification by active cells overtakes signal dissipation in the passive medium. In this section, we examine this process in further detail and derive equations that can be used to determine the frequency and amplitude of auto-induced oscillations when the instability takes place. For simplicity, we shall consider a situation where diffusion of the signalling molecules in the medium is very fast so that spatial variations of s is suppressed. Consequently, the notion of a well-defined transition to the oscillating state can be introduced.

The phase matching condition and threshold cell density
Given that individual cells couple to each other only through the signal field s, a self-consistency procedure similar to the solution of mean-field models in statistical physics can be employed. In this case, the linear equations governing an eigenmode with eigenvalue λ can be divided into subgroups associated with individual cells. The internal variables of a given cell appear in one and only one of the subgroups. Solution of the subset of equations for cell j yields the cell activity ã j =R a,j (iλ) s . The functionR a,j (iλ) is the same function introduced in the preceding section to describe the linear response of a j to a sinusoidal perturbation at frequency ω = iλ. Likewise, a response functionR s (ω) from the linearised relaxational dynamics of s can be obtained, treating contributions from cells as source terms, as in Eq.
(1) of the Main Text. Combining the two steps, we arrive at the following eigenvalue equation for λ, When a particular eigenvalue crosses the imaginary axis, its real part vanishes, while its imaginary part ω o (the onset frequency) satisfies, HereRā(ω) ≡ N −1 N j=1R a,j (ω) is the averaged single-cell response function. Equation (11) can be written separately for the phase shift φ = − argR and amplitude |R| of the response functions. For α 1 > 0, we have, Eq. (12a) determines the frequency ω o at the onset of collective oscillations, while Eq. (12b) gives the threshold cell density N o . As we mentioned in the Main Text, when the signal is passive, phase lead by the cell is required for Eq. (12a) to be fulfilled. The explicit relation presented here complements the energy argument based on Eq. (7), with the activity-generated thermodynamic force O t being proportional to α 1 a j . As it stands, the cell density N does not appear explicitly in the phase-matching condition (12a). Therefore the frequency of collective oscillations can be estimated from separate measurements of the single-cell response and the medium response. In reality, it is conceivable that properties of the medium are affected by the presence of cells, e.g., the concentration of the signalling molecules secreted. Consequently, bothR s (ω) andRā(ω) may have certain weak dependence on N .

The amplitude equations and frequency shift
Beyond the initial instability, nonlinear effects need to be treated explicitly to determine the amplitude and frequency of oscillations. Assuming a periodic state, the signal strength s(t) can be expressed as a Fourier series that includes the first harmonic as well as higher order harmonics produced by nonlinearities in the system dynamics. Likewise, the noise-averaged cellular activity a j (t) can also be expressed as a Fourier series in t with the same basic frequency. For weak noise, the trajectory of the system falls on a well-defined limit cycle whose mean radius r sets the overall amplitude of oscillations, while the amplitude of the nth order harmonic scales as r n . This structure allows for a systematic determination of the amplitudes using perturbation theory. Below, we illustrate the procedure in the case of cubic nonlinearities in both the dynamics for s and the dynamics for a, and comment on similarities and differences in more general situations. When the cell's activity is noisy, more sophisticated schemes based on the probability distribution function of the cellular state need to be introduced (see, e.g. Ref. 11 ).
Let us consider a noiseless version of the adaptive dynamics defined in the Main Text ( Fig. 3 and Methods), together with a modified version of Eq. (1) that includes a cubic nonlinearity, Here τ s = γ s /K s gives the relaxation timescale for the signal. We also set α 1 → K s α 1 for notational simplicity. The two coefficients c 3 and d 3 set the strengths of nonlinearities in the cellular and signal dynamics, respectively. The model has the inversion symmetry s → −s and (a j , y j ) → (−a j , −y j ), all j. Furthermore, if we redefine the sign of s and at the same time change the sign of α 2 and α 1 , the equations remain the same. We now seek a periodic solution to Eqs. (13) in Fourier form, The amplitudes and phase shifts, all assumed to be real, satisfy a set of equations which can be derived by substituting Eqs. (14) into Eqs. (13), and grouping terms according to the order of the harmonic.
Starting from the first harmonic in the expressions (14), the cubic terms in Eqs. (13a) and (13c) generate the first and third order harmonics according to the identity (cos φ) 3 = (3 cos φ + cos 3φ)/4. Hence terms such as A 3 j and B 3 are present in the equations for the first harmonic. On the other hand, the cubic nonlinearities do not generate even order harmonics if they are not included in the series initially. Hence, up to the third order in the amplitudes, the equations for the coefficients of the first harmonic take the form, Here we have introduced the short-hand notationss = B,ã j = A j exp(−iφ a,j ), andỹ j = C j exp(−iφ y,j ).
To gain an intuitive understanding of the oscillatory solution as the cell density increases beyond the threshold N o , we first eliminate the intracellular variableỹ j in Eqs. (15a) and (15b) to obtain, whereR + a,j (ω) ≡ã is a "nonlinear response function" which expresses the ratio of the complex amplitudes of the first harmonic on the limit cycle. Similarly, Eq. (15c) can be rewritten as whereR + s (ω) is a "nonlinear response function" of s on the limit cycle. It is easy to see thatR + a (ω) andR + s (ω) reduce to their respective linear counterpartsR a,j (ω) andR s (ω) when the oscillation amplitudes vanish.
We now combine Eqs. (16) and (18) to obtain the self-consistency condition, which is reminiscent of Eq. (11). HereR + a (ω) ≡ N −1 N j=1R + a,j (ω) is the averaged single-cell nonlinear response function. When all cells are identical,R + a (ω) =R + a (ω). As before, Equation (20) can be rewritten in terms of the phase and amplitude of the nonlinear response functions, Formally, Eq. (21a) can be used to determine the frequency shift at a finite amplitude of oscillation, while Eq. (21b) relates the oscillation amplitude to the cell density N . Since the amplitudes enter quadratically into the nonlinear response functions, they are expected to increase as (N − N o ) 1/2 just above the threshold cell density N o , e.g., the transition is of the Hopf bifurcation type.
In the Main Text, we have considered the case c 3 = 1 and d 3 = 0. Numerically, the oscillation frequency is found to decrease as the coupling strengthN ≡ α 2 α 1 N increases (see also Supplementary Figure 1A). This is consistent with Equation (21a) whose solution at selected oscillation amplitudes is shown in Fig. 3F in the Main Text. As the amplitude of the oscillations increase, φ + a (ω) decreases on the low frequency side. Consequently, the intersection point with φ + s (ω) = φ s (ω) shifts to lower frequencies. Interestingly, the limit cycle associated with the oscillating state in this model shrinks to a fixed point whenN exceeds an upper threshold valueN b . The dependence ofN b on the adaptation error , which is assumed to be small, can be estimated as follows. At a fixed point of Eqs. (13) at d 3 = 0, we have s α 1 N a (fromṡ = 0), y α 2 s (froṁ a = 0), and a y (fromẏ = 0). Consequently, the upper threshold for oscillations has the scalinḡ Figure 1B shows the numerical values forN o andN b against the adaptation error obtained in our simulations, which confirms (22). The oscillating state expands over a larger range of cell densities when individual cells are more adaptive.
Next, consider the case of nonlinear signal relaxation (d 3 = 1) and linear adaptation (c 3 = 0). The onset of collective oscillations is similar to the previous case ( Supplementary Figure 2A), except that oscillations speed up as the cell density increases further. From Equation (19), we obtain which decreases as the oscillation amplitude increases. As shown in Supplementary Figure 2C, the intersection point shifts to the right. The predicted signal oscillation amplitude B = |s| and frequency shift agree quantitatively with our numerical results (Supplementary Figure 2C).
In the more general case when both c 3 and d 3 are nonzero, we need to first expresss andã j in terms of a common variable that specifies oscillation amplitude before applying the phase-matching condition Equation (21a). As we see from the discussions above, depending on which of the two cubic nonlinearities is stronger, the oscillation frequency may shift either to lower or higher values. In general, nonlinearities may also be present in the dynamics of other intracellular variables which need to be dealt with case by case.
When quadratic nonlinearities are present in the system dynamics, the second harmonic is generated and need to be considered in the perturbative analysis. Consider for example the equation for a j with an extra term c 2 a 2 j . Following the same procedure that led to Eqs. (15), we find an additional term c 2ã * jã j on the right-hand side of Equation (15a), whereã (2) j is the amplitude of the second harmonic in a j (t) (including phase). The amplitude equation for the second harmonic relatesã (2) j to c 2ã 2 j and α 2s (2) . Together with the equation fors (2) , amplitudes of the second harmonic can be expressed as a linear combination of terms c 2ã 2 j from different cells. The upshot of this exercise is that coefficient of the cubic term |ã j | 2ã j in Equation (15a) should contain additional contributions proportional to c 2 2 . The nonlinear response functions (17) and (19) on the limit cycle can still be defined in the same way, and Equation (20) still holds formally. Throughs (2) , terms |ã k | 2 from other cells enter the expression forR + a,j (ω). Two conclusions can be drawn from this fact: i) as in the case of cubic nonlinearities, the transition is still of the Hopf bifurcation type; ii) R + a,j (ω) can no longer be determined by simply measuring the response of a given cell to a sinusoidal stimulus at finite strength, as it is affected by the oscillation pattern of other cells in the system due to the quadratic nonlinearity.

Supplementary Note 3. Minimal rate of signal decay/clearance at the onset of collective oscillations
Experiments have indicated that sufficiently fast breakdown of the signalling molecule is needed for DQS in dicty 12 and for sustained oscillations in yeast cell suspensions 13 . Below we derive an upper limit for the signal relaxation time τ s to satisfy the phase matching condition Equation (4) in the Main Text. The result is inversely proportional to the adaptation error of the intracellular circuit.
Taking Equation (6) for the signal phase shift, φ s = − tan −1 (ωτ s ), we see that a longer τ s yields a larger signal delay |φ s | at a given frequency. This is illustrated by Supplementary Figure 3A, upper panel, where the horizontal frequency axis is shown on logarithmic scale. For a given , the intersection of the two phase-shift curves moves to the left, yielding a lower onset oscillation frequency ω o and a larger coupling strengthN (Supplementary Figure 3B). When τ s reaches beyond an upper limit τ * s ( ), the solution disappears (Supplementary Figure 3B). Interestingly, a reduction of in Equation (10) increases φ a (ω) on the low frequency side, and rescues the solution (Supplementary Figure 3A, lower panel).
The observed inverse dependence of τ * s on in Supplementary Figure 3C can be justified from the behaviour of the two phase-shift functions at low frequencies. Close to ω = 0, Equation (9) in the Main Text yieldsR a (ω) ∝ whilẽ R a (ω) ∝ ω, asR a must be an odd function of ω. This is confirmed by expanding Equation (10) in the Main Text at ω = 0. Consequently, φ a (ω) −R a /R a ∝ ω/ . On the other hand, φ s (ω) ≈ −ωτ s in this regime. The two curves has an intersection at low frequencies provided Using the explicit expressions Eqs. (3) and (10) in the Main Text, we obtain τ * s τ y / for small . The critical cell density and onset oscillation frequency at this maximal τ s are given approximately by N o K(α 1 α 2 ) −1 and ω o τ −1 y , respectively. This is compared to N o K(α 1 α 2 ) −1 and ω o (τ a τ y ) −1/2 at τ s < (τ a τ y ) 1/2 , which are insensitive to as long as it is sufficiently small.  Figure 5C, blue box). Below, we present response properties of the model using ACE concentration as the second control variable, in addition to the extracellular glucose concentration. Time is measured in minutes and concentrations in mM.
To move out of the oscillatory regime, we lower the mean acetaldehyde concentration to ACE 0 = 0.05. Experimentally, this can be achieved by adding cyanide (KCN) which reacts with ACE in the solution 17 . Supplementary Figure 6 shows the time course of metabolites under a step-wise increase in the ACE concentration. The four panels are organised following the order of metabolites along the glycolytic pathway, with the addition of ATP, ADP and NAD. Most metabolites adapt at least partially, except F16P and TRIO upstream of the reaction GAPDH that uses NAD and NADH as cofactors. The redox pair NAD and NADH, being tightly connected to ACE, do not adapt either.
We now consider oscillations of the same set of metabolites stimulated by a periodic redox signal at the frequency ω 0 of spontaneous oscillations mentioned above. In Supplementary Figure 7A, ATP, G6P and F6P are approximately in phase with each other, but they are out of phase with Glci at the entry point of the pathway. The non-adaptive F16P has a behaviour of its own. The phase relations for these metabolites have been measured experimentally, and the results agree well with our numerics 18 . In Figs. 7B-C, metabolites from BPG down to PEP share nearly the same phase with each other and with ATP. The non-adaptive TRIO lags slightly behind F16P. In Supplementary Figure 7D, PYR at the end of the glycolytic pathway has an approximately 90 • phase lead over ATP, and furthermore a smaller phase lead over ACE and NAD.
Supplementary Figure 8 shows the phase shifts of metabolites against a sinusoidal signal ACE obtained from our numerical simulations over a broad frequency range. Apart from PYR, the phase relationships among metabolites at ω 0 hold also at lower frequencies. In Supplementary Figure 8B, it is seen that NAD has essentially the same phase as ACE in the frequency interval, while NADH is completely out of phase. Therefore, on the timescale τ 0 , the phase information of ACE is passed without delay onto the redox ratio NAD/NADH, and fed into the network through the reaction GAPDH. Around ω 0 , the phase lead of NADH over ACE is slightly below 180 • , as observed in experiments on glycolytic oscillations 19 . Supplementary Figure 8C shows the downstream metabolites from BPG to PEP oscillate in phase with each other for ω ≤ ω 0 , meaning the internal time scales for this part of the pathway are shorter than τ 0 . In contrast, PYR develops a phase lead in the intermediate frequency regime, as indicated by the two black arrows in Supplementary Figure 8C. (Note that in Fig. 5D of the Main Text, the phase lead extends to zero frequency indicating that the width of the regime depends on the glycolytic flux.) The adaptive variable ATP also has a phase lead in the entire low frequency region.
Supplementary Figure 8D shows the following phase relations between ATP and several other metabolites as summarized by the equations below, The first two relations among the nucleotides ATP, ADP and AMP simply reflect the conservation of their total number, and that the fraction of AMP is much lower than the other two. In-phase relations apply to substrates BPG and PEP of the ATP harvesting reactions PGK and PYK, respectively, while the out-of-phase relation is observed for GLCi in the ATP consuming reaction GLK. The fact that these relations hold almost strictly in the entire frequency region suggests that quasi-steady-state conditions apply to these and neighbouring reactions. It also suggests a prominent role of ATP in synchronising the phase of metabolites distributed along the glycolytic pathway.
In summary, our numerical results suggest the following mechanism of adaptation. Under a stepwise increase of ACE concentration, the information is passed with negligible delay to the redox ratio NAD/NADH, and then through the delayed reaction GAPDH to BPG and PEP, transiently boosting ATP production. The transient increase of ATP concentration then reduces the upstream glycolytic flux by inhibiting the reaction PFK, which in turn decreases the downstream TRIO concentration, eventually returning the GAPDH flux to its pre-stimulus level. Although many metabolites adapt, the negative feedback loop of ATP production appears to be the core. Fig. 5B in the Main Text shows a more complete phase diagram of the response properties at other values of ACE 0 and Glco concentrations, including the region of spontaneous oscillations.

Supplementary Note 5. A minimal model for ATP autocatalysis coupled to intracellular redox state
We constructed a minimal model to test various quantitative aspects of the adaptation mechanism described above. Reduction in the number of dynamic variables is achieved by lumping consecutive metabolites along the linear pathway that are phase synchronised into a single variable denoting their total concentration. This is a reasonable approximation when interconversion among these metabolites is much faster than the time of interest, e.g. the oscillation period. Supplementary Figure 9A illustrates the selected variables and their interactions. Here, y represents intermediate metabolites that do not adapt (F16P and TRIO), thereby playing the role of a memory node. The variable z represents metabolites from BPG to PEP along the glycolytic pathway. The ATP concentration is denoted by p, while the concentration of PYR, substrate for the ACE producing reaction PDC and thus the corresponding cell activity here, is denoted by a. Since NAD (NADH) is always in phase (out of phase) with ACE, we adopt the redox ratio NAD/NADH as the signal s instead. Motivated by a phenomenological two-component model for glycolytic oscillations in Ref. 20 , we introduce a minimal model of glycolysis with redox control as follows: Here, 2p/(1 + p 2h ) gives the reaction flux of PFK that consumes ATP and is also inhibited by ATP at high concentrations (i.e., the negative feedback loop), with the inhibition strength set by the exponent h(> 1/2). The entry carbon flux into the glycolysis pathway is assumed not to be rate limiting, e.g.., one is in a situation of high extracellular glucose concentration. The term (α 2 s + c 0 )y gives the reaction flux of GAPDH, where c 0 sets the "basal" enzyme velocity at s = 0. Leakage of TRIO into the side branch is represented by y. The term 2z/(1 + p 2 ) gives the reaction flux of PYK (and also PGK), which produces ATP but is also inhibited by ATP. In Equation (26c), the stoichiometric factors 1 and 2 in the first two terms on the right-hand-side correspond to the ATP consumption and production upstream and downstream of TRIO, respectively. ATP consumption by the cell outside of glycolysis (e.g., ATPase activity) is modelled by the term 2p 2 /(1 + p 2 ), which grows with the ATP concentration until saturation at a maximal value 2. The output variable a (PYR) is produced by the same flux that produces p (ATP) and degraded at a constant rate α 1 , The glyoxylate shunt GLYO, which is active at ACE 0 0.2 mM or above in the full model, is turned off. For simplicity, we have chosen the time constants on the left-hand-side of the equations to be the same. As we show below, this choice is adequate for recovering the main low frequency properties of the full model.
Supplementary Figure 9B shows the response of the dynamical variables to a sinusoidal redox variation centred around s 0 = 0.04. Except the buffer variable y, all other variables show adaptive behaviour, with a gaining a phase lead of 90 • over p. For comparison, we show in Supplementary Figure 9C  We note in passing that the two-component model of Chandra et al. 20 also exhibits spontaneous oscillations when the rate constant k of the pyruvate kinase reaction (PYK in Supplementary Figure 9A) takes on intermediate values.
As k affects the delay time of the negative feedback control in ATP production, in this sense it plays a similar role as s 0 . However, our model contains an additional buffer node TRIO which is necessary for the adaptive behaviour seen in Supplementary Figure 9. We have also made the ATP consumption rate dependent on the ATP concentration to eliminate certain pathological aspects of the Chandra et al. model at low values of p. Furthermore, our numerical analysis suggests that a sufficiently small but finite adaptation error associated with low flux diversion is needed to reproduce the response diagram Supplementary Figure 9D. On the high (oxidative) end of s 0 , the reaction GAPDH drives down y (TRIO) and hence the flux of the side reaction, making the system adaptive even when ∼ 1.
Comparing the response diagrams of the minimal model (Supplementary Figure 9D) and of the full model at high extracellular glucose concentrations (Main Text, Fig. 5B), we see that the adaptive regime on the oxidative side is restricted to a much narrower region in the latter case. Upon a detailed investigation of the full model we found that, at higher values of ACE 0 , the side reaction GLYO is activated. Shutting down the reaction, we obtained a response diagram similar to that of the minimal model (Supplementary Figure 10). The reaction GLYO uses NAD as cofactor and consumes ATP (see Supplementary Figure 4). With regard to the change in NAD/NADH ratio upon an upshift of ACE, it has an opposite effect as compared to ADH. This and ATP consumption by GLYO leads to a sign reversal in the transient response of ATP to ACE upshift at ACE 0 0.2 mM in the full model (Fig. 5B in the Main Text). Inhibition of GLYO eliminates the sign switch and makes PYR adapting to ACE over a much larger region of the phase diagram (Supplementary Figure 10).
Supplementary Figure 11 shows representative time courses of the PYK reaction flux to a stepwise ACE signal, computed using the original and modified glycolysis model, as well as the minimal model. Concentration of its product, PYR, is found to be proportional to the PYK reaction flux in all three models, i.e., the degradation rate of PYR is a constant. The original and modified models exhibit nearly identical adaptive response on the low ACE (reductive) side, but differ on the high ACE (oxidative) side. In the latter case, the PYK flux is significantly higher and also non-adaptive when the glyoxylate shunt (GLYO) is on. Further numerical investigations of the full model with blocked GLYO reaction show that it shares the following features of the minimal model as the oxidation level increases: 1) the frequency inside the oscillatory regime increases; 2) (mean) p (ATP) and z (BPG, P3G, P2G, and PEP concentrations) increase by a moderate amount; 3) y (TRIO and F16P concentrations) decreases; 4) a (PYR concentration) first increases, then decreases. Experimental time-course measurement with blocked glyoxylate shunt will serve to validate or improve the model assumptions.

Supplementary Note 6. A model for glycolytic oscillations in yeast cell suspensions
To study collective oscillations in a population of cells whose internal dynamics follows Eqs. (26), we adopt the following signal dynamics as in Ref. 21 : Here s in and s ex are the intracellular and extracellular signal concentration, respectively; D is the membrane permeability of the signalling molecule; k in and k ex are the intracellular and extracellular signal degradation rate; and ρ is the volume fraction of yeast cells, which increases with the cell density, and saturates at 1. Hereafter, we call ρ the cell density directly. In the steady-state, the extracellular signal strength (i.e., acetaldehyde concentration) is a function of ρ and D.
Let us first consider the situation of fast equilibrium between s in and s out . Previously, Silvia De Monte et al. proposed a diffusion timescale τ s ≈ 0.003 s by assuming a quasi-stationary concentration profile and that ACE molecules need to diffuse across a spherical shell with an inner radius r 1 = 3 µm and an outer radius r 2 = 6.5 µm 22 . This diffusion timescale is much smaller than the oscillation period of 37 s. Assuming the time for an ACE molecule to cross the cell membrane is of the order of 1 s or less, we obtain the following approximate equation for s = (s in +s ex )/2, Up to corrections of order , the stationary state of the dynamical system defined by Eqs. (26) and (28) is given approximately by, The signal strength increases with the cell density and saturates at 1/(k in + k ex ). At small but finite , corrections to the above expressions become significant at large y or small s, where the side reaction G3PDH in Supplementary  Figure 4 is activated to divert the glycolytic flux. In the numerical studies presented below, we set the two ACE degradation rates k in and k ex to be small. The signal strength s varies over a broad range as the cell density ρ increases.
Supplementary Figure 12 shows numerical solutions of the coupled minimal model at four selected ρ values. Except the case at ρ = 0.01, oscillations of s and the intracellular variables are seen. In Supplementary Figure 13, we plot the oscillation amplitudes and time-averaged values of s and a against the cell density ρ. From the lower panel of Supplementary Figure 13A we see that, for the signal dynamics chosen, the lower adaptive regime in Supplementary Figure 9D is mapped to a narrow interval of cell density 0.003 < ρ < 0.013. From Supplementary Figure 12, we see that onset of collective oscillations in the coupled system takes place somewhere between ρ = 0.01 and 0.014. More detailed studies indicate that the transition is not the expected Hopf bifurcation type, but instead emergence of a limit cycle at finite amplitude. Similar behaviour was seen in the study of the full kinetic model (see Fig. 10 in 16 ). On the other hand, experimental work seem to support the Hopf bifurcation scenario 22,23 . The discrepancy could be due to cell-to-cell variability in the experimental system.
Beyond the onset point, oscillation amplitudes vary continuously with the cell density. For ρ > ρ c = 0.34, the time-averaged value of s falls in the upper adaptive regime in Supplementary Figure 9E. Since the cell density here already exceeds the threshold value required for collective behaviour of adaptive units, oscillations continue.
Finally, we present numerical results demonstrating the effect of a slower cross-membrane transport of acetaldehyde on the collective dynamics. The system dynamics is defined by Eqs.  Figure 14 shows the oscillation amplitude of s in together with the time-averaged values of s in and s out at selected values of D. At D = 100, s in and s out are nearly identical and the system behaviour is essentially the same as described above under the fast equilibrium assumption. At D = 1, the time-averaged value of s ext is noticeably smaller than that of s in , indicating a significant gradient of acetaldehyde concentration across the cell membrane. Nevertheless, collective oscillations via DQS continue at all cell densities. As we further reduce D, DQS stops at large cell densities, as illustrated by the case for D = 0.4. Collective oscillations disappear at D = 0.1. Here, s in remains high due to the slow intracellular degradation rate k in , which places the single-cell dynamics in the upper adaptive regime even when the cell is isolated. However, the phase delay across the cell membrane changes the response properties of the cell to external signal variations. In this case, s in should be considered as the sender of the external signal but as one can see from Eq. (27a), the adaptation (phase lead) of a to s in does not translate to phase lead of s in to s ex when D is small. The latter is required for the adaptation route to collective oscillations. Now, we apply our linear response theory to predict the disappearance of DQS in this system at D = 0.4, where the oscillation amplitude decreases continuously to zero at large cell densities (Supplementary Figure 14). The onset condition cannot be predicted from this theory as collective oscillations emerge first through synchronisation of individual oscillators. Mapping this model to our generic framework, the intracellular signal s in is the "activity", while the extracellular signal s ex is the mediating signal. Absorbing the contribution of cell density into the signal response, i.e.,R sex (ρ, ω) ≡s ex s in = ρD k ex − iωτ s + ρD , as computed from Equation (27b) using Fourier transform, the self-consistency for collective oscillations in this example requires φ sin (ρ c , ω c ) = −φ sex (ρ c , ω c ), |R sin (ρ c , ω c )R sex (ρ c , ω c )| = 1.
A simple form of signal relay efficiency no longer exists due to the nonlinear effect of cell density on the signal spectrum.
Here,R sin ≡s in /s ex measures how an individual cell, defined by Equation (26) and Equation (27a), responds to a periodic perturbation of s ex at frequency ω. As glycolysis is a nonlinear process,R sin also depends on the average levels ex of extracellular signal at which small periodic stimulation is superimposed. As we are more interested in the response spectrum at its natural context of interacting populations, we perturb the cell ats ex determined from the population dynamics at a given cell density ρ. Hence, the signal-level dependence ofR sin translates into celldensity dependence. Given the numerically computed two-dimensional spectrumR sin (ρ, ω) andR sex (ρ, ω), the critical (fixed) point ρ c can be found via numerical iteration. At the critical cell density ρ c , we expect the two equations in Equation (31) to be satisfied at the same frequency. This is indeed observed in Supplementary Figure 15A. Note that the phase shift φ sin is positive only in the intermediate frequency regime, unlike the case for typical adaptive variables seen in other examples in this paper. While φ sin has a phase-leading regime induced by the upper-stream adaptive variable a, s in may not necessarily be adaptive to s ex , especially at large D, where s in and s ex would be almost identical at any time. The predicted critical cell density and frequency agree well with results from direct simulation of the population dynamics (Supplementary Figure 15B). In summary, under fast equilibration between intracellular and extracellular acetaldehyde concentrations, the coupled system exhibits collective oscillations over a broad range of cell densities, encompassing the adaptive and oscillatory regimes of a single cell. Onset of collective oscillations at low cell densities exhibit complex behaviour due to the assumed sensitivity of the reaction GAPDH to the NAD/NADH ratio. Delay in the cross-membrane transport of acetaldehyde weakens adaptation of intracellular metabolite concentrations to change in the extracellular acetaldehyde concentration, and may eliminate collective oscillations altogether when the delay is too long 13 . At moderate delays, rise in the intracellular acetaldehyde concentration brings individual cells to the oscillatory regime even when in isolation. The enhanced oscillation amplitude at D = 1 and low cell densities seen in Supplementary Figure 14, however, is obtained under the assumption that all cells in the population behave identically. This behaviour is susceptible to cell-to-cell variations as well as temporal noise in intracellular dynamics. As D decreases further, an isolated cell exits from the oscillatory state and moves into the upper adaptive region, as shown in the last two columns in Supplementary Figure 14. Collective oscillations at D = 0.4 is attributed entirely to DQS. Our model study exposes this and other subtleties that can affect interpretation of collective oscillations. The specific effects we identified in this work could serve to guide the design of future experiments where various model parameters can be controlled quantitatively, e.g., k ex for extracellular degradation rate of acetaldehyde by adjusting the flow rate in microfluidic setups 17 .