The Interplay between Feedback and Buffering in Cellular Homeostasis

Buffering, the use of reservoirs of molecules to maintain concentrations of key molecular species, and negative feedback are the primary known mechanisms for robust homeostatic regulation. To our knowledge, however, the fundamental principles behind their combined effect have not been eluci-dated. Here, we study the interplay between buffering and negative feedback in the context of cellular homeostasis. We show that negative feedback counteracts slow-changing disturbances, whereas buffering counteracts fast-changing disturbances. Furthermore, feedback and buffering have limitations that create trade-offs for regulation: instability in the case of feedback and molecular noise in the case of buff-ering.However, becausebufferingstabilizesfeedback and feedback attenuates noise from slower-acting buffering, their combined effect on homeostasis can be synergistic. These effects can be explained within a traditional control theory framework and are consistent with experimental observations of both ATP homeostasis and pH regulation in vivo . These principles are critical for studying robustness and homeostasis in biology and biotechnology.


In Brief
Hancock et al. demonstrate that buffering, the use of reservoirs of molecules to maintain stable concentrations of key molecular species such as ATP, and negative feedback can work hand-in-hand to carry out robust cellular homeostasis.

INTRODUCTION
Homeostasis and robustness are fundamental to life (Cannon, 1932;Kitano, 2004Kitano, , 2007Sherwood, 2013). Without them, living organisms can quickly become diseased or die. For a cell to maintain normal function, its levels of pH, ATP, metabolites, and ions, among other factors, must be kept within a tight range, despite continual environmental perturbations and variability in biochemical processes. Studying the principles that enable such robust regulation is critical to topics ranging from the progression of complex diseases and the design of new therapeutic treatments (Kitano, 2004)  Negative feedback and buffering are the primary mechanisms behind the robust regulation of molecular concentrations in biochemical processes (Kitano, 2004;Alon, 2007;Iglesias and Ingalls, 2009;Ma et al., 2009;Sherwood, 2013). In negative feedback, the concentration of a ''regulated'' molecule is sensed, and the information is used by a biological actuator in such a way as to ensure that any concentration changes of the regulated molecule are opposed (Sherwood, 2013). Biological actuators come in many shapes and sizes, from a small molecule effector changing the activity of an enzyme, to the respiratory muscles carrying out breathing. For example, negative feedback counteracts a decrease in cellular ATP by increasing the activities of metabolic enzymes in glycolysis.
Buffering provides an alternative regulatory strategy to negative feedback. A buffer can be defined generally as a reservoir of molecules used to maintain the concentrations of regulated molecular species. In this paper, we specifically consider a reservoir of buffering molecules connected to a pool of regulated molecules via non-actuated, Le Chatelier-driven interconversion. For example, a weak acid (HA) acts as a reservoir to regulate pH (H + ), achieved via acid dissociation (HA4H + +A À ). There is overlap between the actions of feedback and buffering. For example, both Le Chatelier's principle and negative feedback achieve robust regulation by counteracting changes. However, in this paper, we consider cases where the two are separated biologically, based on whether the regulation is via actuation or not.
Despite the fact that feedback and buffering are both often found to be regulating the same processes simultaneously, to the best of our knowledge, there have been no studies investigating the universal principles behind how they function together. For example, ATP concentration is regulated in vertebrates by creatine phosphate buffering and by feedback on metabolic pathways such as glycolysis (Berg et al., 2007;Baker et al., 2010). Furthermore, the two mechanisms often appear to act in concert-creatine phosphate initially supplying muscles with energy during strenuous exercise and feedback controlling the energy supply from metabolic pathways over the longer term (Berg et al., 2007;Baker et al., 2010).
In this paper, we use minimal-model analysis to determine the cooperation and trade-offs between negative feedback regulation and buffering regulation in biochemical systems. We find that feedback and buffering work hand-in-hand to regulate against disturbances (perturbations or unwanted signals) of different frequencies. We also find that buffering acts to improve stability-imposed limits on feedback regulation, but that rapidly acting buffering can introduce high-frequency molecular noise. When buffering is slower acting, however, feedback can attenuate this noise, creating a mutualistic picture where feedback aids buffering with noise and buffering aids feedback with stability. Furthermore, we show the equivalence of rapidly acting buffering with negative derivative feedback control, which is widely used in technology. We also show that slower-acting buffers can act as a barrier or filter between the regulated species and disturbances in the buffering species. We take our general findings and apply them to the specific case of cellular ATP regulation, in order to provide new insight into the cellular strategies used for energy management.
The use of a minimal model, independent of mechanistic details, provides simplicity and allows us to more easily highlight the underlying principles behind regulatory behaviors. Furthermore, our definition of buffers-reservoirs of molecules used to maintain concentrations of other molecules-encapsulates a wide array of reservoir arrangements, from spatial isolation of molecules by membranes and compartments to the sequestration of macromolecules by non-covalent binding. In addition, the buffer systems that we consider are open systems that continuously exchange mass with their environment. The familiar closed-system treatment of buffers (e.g., pH buffering, where the disturbance to the system is a one-time addition of a chemical) represents a limit of the more general open-system problem studied here.

RESULTS
We begin by introducing a minimal model that consists of both a regulated and a buffering biochemical species. Illustrated graphically in Figure 1A, the model can be represented mathematically by _ y = p y ðyÞ |fflffl{zfflffl} production with feedback À g y ðy; xÞ + g x ðx; yÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} buffering À g y y |{z} removal + d y |{z} disturbance _ x = g y ðy; xÞ À g x ðx; yÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} buffering ; (Equation 1) where y is the regulated species, x is the buffering species (the ''reservoir''), p y is the production rate of y, g y is the kinetic rate constant for the removal of y, the buffering reactions are g y (y,x) representing y/x and g x (x,y) representing x/y, and d y is the disturbance. Incorporation of feedback into the system is represented by the y dependence of the production term p y ; feedback senses the regulated variable y and actuates the production rate p y to compensate for changes. For simplicity of exposition, we restrict our focus to feedback exclusively on the production of species y. This model can be used to represent numerous biochemical systems. For example, it can serve as a simple model of ATP regulation ( Figure 1B), where y is ATP concentration, x is creatine phosphate (pCr) concentration, the production p y represents glycolysis, the removal g y represents nominal ATP demand, and d y can be used to represent a range of disturbances, from upstream nutrient disturbances to changes in ATP demand in muscles during exercise.
We define the ''nominal steady state'' of the system to be its zero-disturbance (d y = 0) steady state, and denote it as ðy; xÞ; here y is called the set point of the regulated variable (see STAR Methods A.1). To analyze the regulation properties of the model, we linearize Equation 1 and reformulate it in terms of deviations from the nominal steady state (Dy,Dx) (see Equation 12 in the Appendix and STAR Methods A.2). More general cases that include non-linear removal rates and the direct removal of the buffering species are treated in the STAR Methods and provide qualitatively similar results.
Throughout this work, we have studied both the cases of rapid and non-rapid buffering, where rapid buffering reactions occur on a much faster timescale than all other reactions and so rapidly reach a (quasi) equilibrium (see STAR Methods A.3). In each of the following sections, we initially present results for the rapid buffering limit, as they are mathematically simpler. We then discuss results for the more general non-rapid buffering case, for which the mathematics can be found in STAR Methods A and B.

Feedback Rejects Low-Frequency Disturbances and Buffering Rejects High-Frequency Disturbances
Our first focus is on the forced response of the minimal model to simple oscillating disturbances. This is because any disturbance signal can be decomposed into its fundamental oscillatory components (through the Fourier transform) (Nise, 2004). A constant disturbance has a frequency of zero and so is also covered by this analysis. . y represents the regulated species, x the buffering species, and d y the disturbance. x is the buffer reservoir that compensates for changes to y via Le Chatelier-driven y4x conversion. Feedback, represented by the dashed line, senses the y amount and actuates the y production rate, p y , to compensate for changes.
(B) Minimal model for ATP regulation. ATP is the regulated species, pCr is the buffering species. ATP is supplied via the glycolysis pathway and consumed via downstream ATP demand. Perturbations to the supply and demand of ATP represent disturbances to the model.
To quantify the effect of disturbances on the regulated species, we introduce the sensitivity function (see the Appendix and STAR Methods A.4 and A.5): f u = kDyk pow y 0 kd y k pow p y = 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ð1 + HÞ 2 + ð1 + BÞ 2 u 2 n q ; (Equation 2) where kDyk pow measures the time average of the magnitude of deviation Dy from the set point y (kDyk 2 pow represents the power of the signal Dy; see Equation 14 in the Appendix). Similarly, kd y k pow measures the time average of the disturbance magnitude. Also, p y is the nominal production rate, u n is the normalized frequency of the oscillating disturbance, B is the buffering equilibrium ratio, and H is the normalized feedback gain (see the Appendix). A smaller f u sensitivity corresponds to better disturbance rejection.
From Equation 2, we can see that the contribution of buffering (B) to regulation becomes more important as the frequency of the disturbance increases, whereas for lower-frequency disturbances, regulation becomes dominated by negative feedback (H).
These points are further illustrated in Figures 2A-2F and are also consistent with known biological observations of ATP (Berg et al., 2007) and pH regulation (Sherwood, 2013), where feedback (on metabolism or breathing rate) is observed to regulate slow-changing disturbances (i.e., disturbances composed primarily of low-frequency components) while buffers are observed to regulate fast-changing disturbances (composed primarily of high-frequency components). This is often described in the literature as buffering acting initially and feedback acting over the longer term (Sherwood, 2013).
These points also parallel the high-frequency/low-frequency split in the regulation of industrial chemical processes by feedback and buffering tanks (Faanes and Skogestad, 2003).
Feedback and Buffering Work Together to Reject Impulse Disturbances We next focus on the transient response of the minimal model by examining the effect of an impulse disturbance on the regulated species. The impulse signal is a useful test signal for analyzing the rejection of generic disturbances. This is due to it having a uniform frequency spectrum, i.e., the impulse signal contains all frequencies in equal amounts (Nise, 2004).
The sensitivity of the model in Equation 1 to an impulse disturbance is given by (see the Appendix and STAR Methods A.6) where kDyk L2 measures the total accumulation of Dy over time, b d y is the magnitude of the impulse disturbance at d y , and T u = 1/g y is the timescale of the unregulated system. As before, a smaller f imp sensitivity equates to better regulation, and it is clear from Equation 3 that increasing either or both feedback and buffering improves rejection of impulse disturbances. Simulations for the rejection of an impulse disturbance are shown in Figure 2G. We see that feedback shortens the time to Here, y is the set point of the regulated species, and smaller deviations from y imply better regulation. Note that in the oscillatory cases, temporal plots show the response some time after the onset of the disturbance, once oscillations have settled to a steady state.
(A-F) Low-frequency (A and D), medium-frequency (B and E), and high-frequency (C and F) oscillatory disturbances. For the temporal plots, the disturbance amplitude is 50% of the nominal production rate p y . Feedback is the dominant regulatory mechanism at lower frequencies, while buffering dominates at higher frequencies. Note that as their frequency increases, disturbances become increasingly rejected independent of feedback and buffering (i.e., for all regulation cases). This is due to the damping effect of the y removal term (Àg y y in the model).
(G and H) Impulse disturbance. For the temporal plot, the disturbance occurs at t = 0 with a magnitude of 50% of the set point y. In comparison with feedback, buffering initially works much quicker in moving y back toward its set point but suffers from a longer-lived residual deviation.
(I) Comparing impulse disturbance rejection between rapid and slower buffering. Rejection is improved for faster buffering, especially in the short term.
remove excess y, and so reduces the cumulative Dy. In contrast, buffering very quickly reduces a large initial excess of y by rapid conversion to the buffering species x. However, additional time is required to clear the remaining excess, as the x/y conversion from the enlarged buffer pool compensates for any subsequent removal of y.
If we no longer assume that the buffering reactions are rapid, we find that increasing the speed of the buffering reactions (increasing R in Equation 12 in the Appendix) improves the effectiveness of disturbance rejection in the transient phase (see Figure 2I and STAR Methods A.7). In particular, we find that the rapid buffering special case is optimal for impulse disturbance rejection when the disturbance is acting on y (see STAR Methods A.8).
Buffering Improves Stability-Imposed Limits on Feedback Effectiveness Often, regulatory inputs to a system are subject to delays that limit the effectiveness of feedback regulation. For example, feedback delays occur in blood pH regulation, where the delay is caused by the time taken for blood to be transported around the body. In contrast, buffering mechanisms typically do not contain significant delays. In this section, we show that this property enables buffering to effectively stabilize feedback loops made unstable by delays, thereby improving fundamental limits on the effectiveness of feedback regulation.
To investigate this effect, we modify our model in Equation 1 by introducing a time delay t to the feedback component of the production term p y (y). This approximates a multitude of delay forms, from a single slow process to a long chain of reactions. The delay does not change the steady state of the system, but it does change its dynamic behavior and can cause the system to become unstable and oscillate. In order for the system with delay to be stable, the feedback gain, H, must be constrained such that (see STAR Methods A.9) H < H crit and therefore f u;low > 1 1 + H crit ; (Equation 4) where The increase of H crit due to an increase in B demonstrates that buffers have a stabilizing effect that can allow for stronger feedback and thus better rejection of slow disturbances (see Figure 3). As H nears H crit , the system begins to demonstrate a resonance effect, becoming highly sensitive to oscillating disturbances near a resonant frequency (see STAR Methods A.8). When the critical gain threshold is crossed (H > H crit ), feedback works too fast relative to the duration of the delay, overcorrecting for deviations to the regulated species and causing unstable oscillations (see Figure 3A).
We also find that the stability constraints on feedback strength are further improved when the buffering reactions slow down (decreasing R in Equation 12; see STAR Methods A.9). This creates a trade-off with the effectiveness of impulse disturbance rejection, for which rapid buffering can be shown to be optimal (see STAR Methods A.8).

Rapid Buffering Introduces High-Frequency Molecular Noise
Biochemical reactions are inherently random molecular processes that introduce molecular noise into biochemical systems (Raser and O'Shea, 2005;Klipp et al., 2009;Lestas et al., 2010;Oyarzú n et al., 2015). Negative feedback can reduce molecular noise (Thattai and van Oudenaarden, 2001;Dublanche et al., 2006;Lestas et al., 2010), while non-rapid buffering alone (without feedback) may increase or decrease molecular noise depending upon the correlations of births and deaths of molecular species (Rooman et al., 2014). In what follows, we show that rapid buffering introduces high-frequency noise into a biochemical system with negative feedback.
To study noise, the minimal model in Equation 1 can be modified to capture the discrete nature of molecule numbers and the probabilistic nature of molecular reactions (see STAR Methods A.10). This model provides probabilistic predictions of the number of molecules changing over time. We use a linear noise approximation and assume rapid, linear buffering to derive a function for noise intensity (see STAR Methods A.10). This measures the sensitivity of our minimal model to molecular noise as A B C Figure 3. Buffering Improves Feedback Stability (A and B) Responses to a step disturbance at t = 0 for the minimal model with feedback delays under various buffering-feedback combinations. The disturbance amplitude is 50% of the nominal production rate p y . In (A) (weak buffering), increasing the feedback gain causes oscillations, while in (B) (strong buffering), increasing the feedback gain instead improves steady-state regulation beyond what is possible when buffering is absent. (C) The feedback gain and buffering equilibrium ratio for the six cases in plots (A and B). The shaded area corresponds to the unstable region, where the stability boundary is determined by Equation 4.
where y n represents the number of molecules of the regulated species (see Appendix), and s y n and hy n i are the SD and average of y n , respectively. The average, hy n i, is used as the set point for the stochastic system. Note that here, we analyze noise under an invariant set point for consistent comparison between different (H,B) arrangements. (See STAR Methods A.10 for further details on treating hy n i as invariant, as well as treatment of the nonrapid, non-linear buffering case.) If there is no feedback, then F = 1 and the molecular noise has a Poissonian distribution characterized by s 2 yn = hy n i. This represents the upper noise limit, as F % 1 for all cases of regulation, i.e., for all choices of H or B. From the form of F, we see that increased feedback (H) always reduces molecular noise. On the other hand, an increase in B always increases molecular noise (except in the complete absence of feedback, where F = 1 always and molecular noise is always at its upper Poissonian limit). These effects are illustrated in Figure 4.
In Equation 5, F is broken into two terms: The first term represents low-frequency noise resulting from the production and removal reactions. The second term represents noise created by the buffering reaction, which is high-frequency for rapid buffering (see the Appendix for noise bandwidths). This second term is independent of feedback because the noise frequency is too high to be effectively regulated by feedback. These different noise-frequency characteristics can be observed in the temporal plots in Figure 4.
As the buffering reaction is slowed down, the frequency composition of the buffer-generated noise is lowered (see the Appendix). When the buffering reactions become sufficiently slow, the buffer-generated noise can be effectively regulated by feedback (see STAR Methods A.10). In this case, a synergistic effect emerges, where feedback aids buffering with noise, just as buffering aids feedback with stability.
Rapid Buffering Is Equivalent to Negative Derivative Feedback PID (proportional-integral-derivative) feedback controllers are crucial regulatory components of many technological systems. They provide a straightforward method for considering and combining past (integral), present (proportional), and future (derivative) state information (Aströ m and Murray, 2008). While proportional and integral feedback have been well studied in biological contexts (Othmer and Schaap, 1998;Saunders et al., 1998;El-Samad et al., 2002;Mello and Tu, 2003;Alon, 2007;Ma et al., 2009;Muzzey et al., 2009;Ni et al., 2009;Ang and McMillen, 2013;Briat et al., 2016), significantly less is understood about the role of derivative feedback.
Here, we show that rapid buffering is mathematically equivalent to negative derivative feedback. This finding generalizes a previous in silico observation that creatine acts as derivative control for ATP (Cloutier and Wellstead, 2010 with the buffer equilibrium constant B equivalent to the derivative feedback gain (see also STAR Methods A.11). A mechanistic interpretation of how the buffering reaction implements derivative feedback is included in STAR Methods A.11. There is a key difference between buffer-feedback regulation and standard PID controllers as a result of feedback delays. In Equation 7, the proportional term due to production is subject to feedback delays, while the derivative term due to buffering is not. Delays fundamentally limit the effectiveness of feedback, including the derivative components of standard PID controllers These simulations demonstrate that rapid buffering introduces highfrequency noise, and therefore effective noise reduction is better achieved through feedback devoid of buffering. (Freudenburg and Looze, 1987;Skogestad and Postlethwaite, 2005), whereas this is not the case for buffers (as discussed above). If a hypothetical derivative action occurred through feedback on production, it would also be subject to delays. Thus, rapid buffering can provide a significant advantage over an equivalent derivative feedback mechanism in the presence of delays.
Another key difference between buffering and standard derivative feedback exists in terms of their regulation performance in the presence of noise. The molecular noise of buffering is similar to measurement noise on a derivative controller (see STAR Methods A.11), i.e., noise that directly affects the controller's ability to sense the regulated variable (Skogestad and Postlethwaite, 2005). However, in systems with rapid buffering, the amplification of noise due to an increased buffer equilibrium, B, eventually saturates (see Equation 5). This is in contrast to technological controllers, which have no corresponding limit on their tendency to amplify measurement noise. For this reason, technological controllers based on derivative control need to filter out high-frequency noise components before applying derivative action. The combination of such pre-filtering with derivative control is called ''derivative filtering'' (Aströ m and Murray, 2008). Interestingly, non-rapid buffering reactions are equivalent to derivative filtering as they can be shown to also filter high-frequency components of noise, which for buffered systems is important for low set points hy n i (see STAR Methods A.11 and also Table 1 for trade-offs).

Slow Buffering Filters out High-Frequency Reservoir Disturbances
Disturbances may also act directly on the buffering reservoir (see Figure 5A and Equation 15 in the Appendix). In this case, nonrapid buffering acts as a low-pass filter, preventing high-frequency components of the disturbance, d x , from reaching the regulated species, y. For oscillating reservoir disturbances, the sensitivity function is (see STAR Methods A.12) where u n is the normalized frequency of the oscillating disturbance, R is the normalized buffer reaction speed (see Equation 12), and f u is the original sensitivity function for a disturbance d y , as defined in Equation 2. Thus, when u n (R, regulation properties for oscillatory disturbances at d x and d y are the same. However, as the factor u n /R increases (through a higher disturbance frequency or slower buffering), disturbances at d x experience further filtering (f dx u <f u ). Simulations shown in Figures 5B and 5C illustrate this effect. Figure 5D shows how such filtering affects the response to impulse reservoir disturbances ( Figure 5D).
This manner of d x filtering is a feedforward rather than feedback property of the system; it takes place independently of changes to y and of the reaction g y , with the reaction g x being directly responsible (see STAR Methods A.12). The combination of feedforward filtering with derivative feedback resembles a two-degree-of-freedom controller in technology (Skogestad and Postlethwaite, 2005).
Biologically, the reservoir ''stores'' the disturbance d x , releasing it slowly into the regulated species. The action is similar to a river dam that collects upstream water flow while releasing a much more constant flow downstream. This effect is commonly replicated in industrial chemical processes by the use of large, upstream buffering tanks (Faanes and Skogestad, 2003).

Slow Buffering with Reservoir Feedback Acts as Integral Feedback
It is also possible for the production rate of y to be dependent on feedback from the buffering species (see Figure 5E). This can be modeled by replacing p y (y) with p y (x) in Equation 1.
When the buffering reaction is rapid, this remains equivalent to feedback emanating directly from y, the only difference being a change to the feedback gain (see STAR Methods A.13). However, for non-rapid buffering, changes to x accumulate; the present-value x becomes dependent on past values of y, and the feedback acquires a so-called memory. As buffering slows, such feedback begins to approximate integral feedback (see STAR Methods A.13), improving the rejection of constant and slow-changing disturbances (Nise, 2004).

Example: Creatine Phosphate Buffering Improves Feedback Effectiveness in Glycolysis
As a case study, we consider the buffering of ATP by creatine phosphate (pCr). Phosphate group donation from a pCr reservoir (pCr + ADP4Cr + ATP) allows for the quick regeneration of ATP when there is a sudden increase in ATP demand. Our results show that pCr buffering also enables significant improvement in glycolytic response and the rejection of slow-changing ATP disturbances. Key to this result is the fact that glycolysis is susceptible to oscillations (Keener and Sneyd, 1998; Yang et al., 2008) and is  [ high-frequency molecular noise (up to the rapid buffering limit) ( Figure 4) Y Buffering reaction speed [ rejection of high-frequency and impulse disturbances at d x [ stability (when speed is not too low) Y molecular noise frequency (such that noise can be better regulated by feedback) Y rejection of impulse disturbances at d y ( Figure 2H) subject to a fundamental, stability-imposed limitation on feedback regulation (Chandra et al., 2011), similar to that discussed earlier in this paper. In glycolysis, instability arises primarily due to its autocatalytic nature: where an initial ''investment'' of ATP is required before more ATP can be generated. The delay between investment and return stages can drive glycolytic oscillations via negative feedback, especially under higher feedback gains, which demand larger initial ATP investments (Chandra et al., 2011). Here, we show that pCr buffering provides stability by supplying ATP for the initial investment stage of glycolysis. This relaxes the fundamental feedback limit and enables improved regulation.
To study this, we modify the minimal model from Chandra et al. (2011) to include pCr buffering. As in Chandra et al. (2011), we assume that the total concentration of adenosine phosphates in the cell, including ADP and AMP, remains constant, and that the activating effects of AMP can be modeled as ATP inhibition. By further assuming total creatine to be constant, we simplify the buffering reaction as pCr4ATP. The resulting model is illustrated graphically in Figure 5F and can be represented mathematically by (see STAR Methods B) _ m = fðyÞ |ffl{zffl} PFK À wðyÞm |fflfflffl ffl{zfflfflffl ffl} PK _ y = À q fðyÞ |ffl{zffl} PFK + ðq + 1Þ wðyÞm |fflfflffl ffl{zfflfflffl ffl} PK À g y ðy; xÞ + g x ðx; yÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} where y represents ATP concentration, x represents pCr concentration, and m is a lumped variable representing the concentration of intermediate metabolites. Also, g y (y,x) and g x (x,y) represent reaction rates for ATP4pCr, f(y) represents the flux through phosphofructokinase (PFK), w(y)m represents the flux through pyruvate kinase (PK), V y represents the nominal ATP demand, and q is the stoichiometric investment/return ratio for the autocatalysis of ATP. As in the general analysis of the minimal model in Equation 1, we begin by linearizing Equation 9 and initially assuming rapid buffering. The feedback gain, H PFK , in the linearized model corresponds to the (negative) normalized slope of ATP to PFK activity (see Equation 16 in the Appendix). Stability imposes that the feedback gain and therefore the steady-state sensitivity of ATP to disturbances in its demand are constrained, such that (see STAR Methods B.1) with B representing the buffering equilibrium ratio, H PK representing the (normalized) negative slope of ATP to PK activity, and k representing the ratio of ATP to intermediate metabolites.
The interested reader is referred to STAR Methods B.3 for the analysis leading to similar results when buffering is not assumed to be rapid (occurring for lower concentrations or different isoforms of creatine kinase). Again, we observe that feedback gain and therefore disturbance rejection are stability constrained, but that the feedback (B and C) A slower buffer filters out high-frequency reservoir disturbances, d x , from the pool of regulated species, y. In these plots, the disturbance amplitude is 50% of the nominal production rate p y . (D) A slower buffer causes an impulse reservoir disturbance to be released more slowly into the regulated pool. In this plot, the disturbance occurs at t = 0 with a magnitude of 50% of the set point y.  (2011), we show that our results also hold for models with significantly more complex feedback structures; specifically, for any stable, ''causal'' form of feedback (see STAR Methods B.4). This includes, for example, feedback from ATP to PFK via fructose 2,6-bisphosphate (Cloutier and Wellstead, 2010). This approach ensures that the results depend only upon the properties of buffering, feedback, and autocatalysis, and do not arise from other case-specific details or simplifications within the model.

DISCUSSION
Buffering regulation can be divided into two main strategies. The first uses a reservoir of molecules to compensate, via Le Chatelier's principle, for changes to the regulated species (see Figure 1); e.g., pCr compensating for varying ATP levels. This strategy typically requires a rapid buffering response for optimal regulation. The second strategy uses the reservoir to block upstream or downstream disturbances from the regulated spe-cies (see Figure 5A). This barrier strategy typically benefits from slower buffering speeds. While blocking a disturbance outright may be preferred to compensating for one after it has taken effect, direct disturbances to the regulated variable cannot always be avoided; e.g., ATP demand acting directly on the cytosolic ATP pool. Therefore, both buffering strategies are complementary.
Motivated by the findings from our theoretical analysis of glycolysis, we searched for an empirical relationship between glycolytic feedback via PFK and ATP buffering by pCr (see STAR Methods B). We observed a correlation across multiple species and cell types between the normalized ratio of the change in PFK activity to the change in ATP concentration and the ratio of pCr to ATP concentrations, as shown in Figure 6. These ratios can be interpreted as proxies for feedback gain, H PFK , and buffering ratio, B, respectively (see Equation 16 in the Appendix and STAR Methods B). The empirical relationship therefore displays qualitative consistency with our generic analysis of glycolysis, which showed that increased ATP buffering enables higher-gain glycolytic feedback.
Further support is offered by published studies that have observed glycolytic oscillations after the removal of pCr: in vitro, in cell-free extracts of rat skeletal muscle (Tornheim and Lowenstein, 1973;Tornheim et al., 1991), and in vivo, in rabbit ventricular myocytes (when aerobic metabolism is also inhibited) (Yang et al., 2008). In contrast, glycolytic oscillations are not observed under typical conditions in species that have native pCr buffering levels and glycolytic feedback gains that are both low (e.g., the single-cell microorganisms in Figure 6).
All organisms also experience ATP and phosphate buffering by inorganic polyphosphates (polyPs), polymers that can contain over 1,000 phosphate residues (Kornberg et al., 1999;Achbergerova and Nahalka, 2011;Rao et al., 2009). In particular, yeast and bacteria can accumulate large reservoirs of phosphate residues in polyP form under certain environmental conditions (Kornberg et al., 1999;Achbergerova and Nahalka, 2011). Phosphate transfer between ATP and polyP differs from pCr in that individual polyP molecules change length, but the total polyP concentration does not change unless a new polyP molecule is formed or an existing one is eliminated. This yields an ATP-polyP buffering equilibrium ratio, B, that is orders of magnitude smaller than for ATP-pCr buffering (see Appendix). In fact, polyP buffering more closely resembles proportional feedback characterized by the reaction rate of ATP-polyP conversion, rather than derivative feedback characterized by its equilibrium ratio, B (see Appendix). This form suggests that polyP buffering is better suited for rejecting slower (lower-frequency) ATP disturbances, while pCr buffering is better suited for rejecting faster (higher-frequency) disturbances and stabilizing oscillations (Aströ m and Murray, 2008). Moreover, while polyP buffering is not constrained by stability (in contrast to proportional feedback in glycolysis), it is constrained by the need to avoid polyP depletion and possibly also excessive accumulation. This moderates the ATP-polyP conversion rate and thereby the potential speed of polyP buffering.
The applicability of our general analysis may extend beyond molecular species that have a dedicated buffer. Reversible sequestration of any form, e.g., molecular complexing not for the primary purposes of establishing a buffer, contributes an as proxies for feedback gain and buffering ratio, respectively (see STAR Methods B), we find empirical cross-species support for a correlation between glycolytic feedback via PFK and ATP buffering by pCr. The shaded area in the figure corresponds to the unstable region, where the stability boundary is roughly and qualitatively estimated for illustrative purposes. This correlation, together with the theoretical results of Equation 10, and previously observed experimental cases of glycolytic oscillations induced by the removal of pCr, suggest that buffering enables higher-gain glycolytic feedback. Red data point, measurements and estimates of in vitro systems; blue data points, in vivo measurements and estimates. additive buffering effect toward the free-form species (i.e., B = S i B i for y4x i ). Thus, several simultaneous forms of sequestration may add up to an effective buffer by sharing the regulatory load. This may lead to a ''hidden'' layer of buffering that should not be ignored. Examples include the heterologous pool of proteins that buffer pH in blood, mRNA transcripts buffering the free ribosome pool ($80% of ribosomes are mRNA-bound in bacterial cells; Forchhammer and Lindahl, 1971), and DNA operator sites buffering the free concentration of a transcription factor (Soltani et al., 2015).
The analysis in this work can be extended to a number of other important topics, including spatial regulation, non-local analysis, time-varying or actuated buffering, alternative buffering topologies, control tracking problems (where the set point may vary with time), combination with feedforward regulation (which anticipates expected signals), trade-offs due to metabolic burden, design in synthetic biology (where buffers have already been used to increase the modularity of synthetic ''parts; '' Mishra et al., 2014), and the study of how numerous other metabolites and ions are regulated within the cell.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: where b y and b x are the forward and reverse buffering rate constants, h y is the feedback gain, and g y is the removal rate constants. The model holds locally about the set point, noting that concentrations will stay near the set point with effective regulation (Skogestad and Postlethwaite, 2005). In the main section, we introduce the following lumped parameters: when the buffering reactions are at (quasi-) steady state and Dy is non-zero.
The feedback gain H can be seen to be the log-log sensitivity of changes in y to changes in production.

Sensitivity Function Definitions
The two sensitivity functions used to evaluate regulation effectiveness (f u in Equation 2 and f imp in Equation 3) are based on the square root of the output power and the L 2 norm, which are defined as respectively, where kDyk pow = kDyk if Dy is at steady state (Skogestad and Postlethwaite, 2005). The output power is useful for persistent disturbances, such as oscillations, while the L 2 norm is useful for time-limited disturbances such as impulses. An impulse of unit magnitude (known as a Dirac delta function) refers to a generalized function satisfying R + a Àa dðtÞdt = 1 for all a>0 (Skogestad and Postlethwaite, 2005), where d(t) = N when t = 0 and d(t) = 0 when t s 0.
In Equation 2, the normalized frequency is u n = u/T u , where u is the frequency of the oscillating disturbance and T u = 1/g y is the timescale of the unregulated system. In Equation 3, the timescale is required to normalize the L 2 norm, as it is integrated over time.

Molecular Noise
The relationship between the molecular number y n and the concentration y is y = y n /U, where U is the system volume. For the linear noise approximation, the set point hy n i = Uhyi = Uy.
We use the term ''stationary'' to indicate a probability distribution that does not change with time. The molecular noise for production and removal has an approximate bandwidth between zero and u low =(g y + h y )/(1 + B) while the molecular noise for buffering has an approximate bandwidth between u low and Disturbance at the Reservoir For studying disturbances at the reservoir, we use the model _ y = p y ðyÞ |fflffl{zfflffl} production with feedback À g y ðy; xÞ + g x ðx; yÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} buffering À g y y |{z} removal _ x = g y ðy; xÞ À g x ðx; yÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} where d x represents the disturbance acting directly on the buffering species x.

Glycolysis Example
The glycolytic feedback gain is (see STAR Methods B) where y is the set point of ATP concentration and f(y) represents the flux through phosphofructokinase (PFK) (see Equation 9).

Polyphosphate Buffering Reaction
For ATP buffering by polyphosphate (polyP), the buffering equilibrium ratio is . As polyP lengthening does not change polyP concentration, D[polyP] is smaller than for an equivalent monomer buffer. The net rate of polyP accumulation can be modeled via Michaelis-Menten kinetics as (Dien and Keasling, 1999) ( Equation 17) where g y and g x are the ATP/polyP and polyP/ATP reaction rates, respectively, E is the polyP kinase concentration, and k 1 , k 2 , K 1 , and K 2 are kinetics constants. This net rate is dependent upon the levels of ATP and ADP, but not the concentration or length of polyP, nor the size of the reservoir of phosphate ions. Assuming conservation of total adenosine phosphates and linearization about a set point, it is possible to rewrite Equation 17 as V polyP = Àb y D[ATP], where b y is the linearized rate constant for ATP/polyP conversion. Thus, V polyP , which represents the contribution to dD[ATP]/dt from polyP buffering, is approximated by proportional feedback.  Figure 6 uses experimental data for: rat muscle (Beis and Newholme, 1975;Bosca et al., 1985;Kushmerick et al., 1992), rabbit muscle (Kemp, 1970;Storey and Hochachka, 1974;Pettigrew and Friedens, 1978), rabbit liver (Ennor and Rosenburg, 1951;Kemp, 1970), desert locust (Walker and Bailey, 1969;Beis and Newholme, 1975), lobster muscle (Beis and Newholme, 1975;Sugden and Newholme, 1975), yeast (Kopperschlager et al., 1968), Escherichia coli (Zheng and Kemp, 1992), and Setaria cervi (Sharma, 2011). The normalized ÀDPFK activity/D[ATP] ratio is determined graphically from the slope of ATP concentration versus PFK activity plots when PFK activity is at half maximum (see STAR Methods B.2); the ratio of ATP to pCr concentrations is determined at steady state.

SUPPLEMENTAL INFORMATION
Supplemental Information includes seven figures and can be found with this article online at https://doi.org/10.1016/j.cels.2017.09.013.  where y is the regulated species, x is the buffering species, p y is the production rate of y, n y (y) and n x (x) are the removal rates of y and x, the (lumped) buffering reactions are g y (y,x) representing y/x and g x (x,y) representing x/y, and d y is a disturbance. Incorporation of feedback into this class of systems is represented by the y dependence of the production term p y . We initially restrict our focus to feedback exclusively on the production of species y.

Conceptualization
For the case presented in the main text, we assume linear removal of species y, such that n y (y) = g y y, and no removal of species x, such that n x (x) = 0. In assuming n x (x) = 0, we look at the class of cellular biochemical systems where the time scale is much faster than the growth rate of the cell. This includes many metabolic processes.

A.1. Nominal Steady State
The following section establishes the steady state (set point) of the minimal model. This is required later for the analysis of deviations about the steady state. The section also defines the steady-state production, the steady-state buffer equilibrium ratio, and the steadystate removal constant. These quantities prove useful for non-dimensionalization in later sections.
The nominal steady state ðy; xÞ of Equation (18)  g y ðy; xÞ = g x ðx; yÞ + n x ðxÞ ðbuffering rates at steady stateÞ p y ðyÞ = n y ðyÞ + n x ðxÞ ðproduction=removal rates at steady stateÞ: The solution ðy; xÞ is implicit, and if required can be determined numerically when an analytical solution cannot be found. In the following, we assume that there is one isolated nominal steady state of interest, and study the properties of the solution as it deviates from this steady state.
For later analysis, it is useful to define the steady-state buffering ratio: the nominal production rate with feedback: p y = p y ðyÞ; and the steady-state removal constant (see also Figure S1): An important case, used in later sections, occurs where removal is linear (i.e., n y ðyÞ = g y y and n x ðxÞ = g x x = g x By), leading to G = g y + Bg x :

A.2. Minimal Model Linearization
The following section determines the linearized minimal model-a useful approximation that permits an analytical and more thorough mathematical analysis. We next linearize the non-linear model. To carry this out, we define deviations from the nominal steady state as Linearizing our full model around the nominal steady state, we obtain where h y = À v vy p y ðyÞj y = y ; b y = v vy À g y ðy; xÞ À g x ðx; yÞ Á j ðy;xÞ = ðy;xÞ ; b x = v vx À g x ðx; yÞ À g y ðy; xÞ Á j ðy;xÞ = ðy;xÞ ;

A.3. Rapid Buffering
The following section analyzes the effect of buffering reactions that occur on a time scale significantly faster than the feedback process. The buffering reaction can be considered to rapidly obtain an equilibrium state, which it maintains while the rest of the system changes more slowly. This 'rapid-equilibrium' approximation reduces the model from two variables to one. Conditions for the rapidequilibrium approximation to hold are then obtained using singular perturbation theory. We initially study the special case of rapid buffering, where we assume that the buffering reactions reach their equilibrium much faster than the production reaches equilibrium with removal. We find a simplified model for the rapid buffering case and find two assumptions for rapid buffering to hold. The first 'weak' assumption (see Equation (22)) is used when the disturbance is bounded. The second 'strong' assumption, presented in a later section (see Equation (34)), is used for unbounded impulse disturbances.
We can use the total quasi steady-state approach (Borghans et al., 1996) to determine a reduced model for rapid buffering, where interconversion reactions between species are 'fast' while the production and removal reactions are 'slow'. This approach allows more accurate reduced models with more general assumptions (Pedersen et al., 2007) for enzyme kinetics (Borghans et al., 1996), signalling network (Pedersen et al., 2007) and gene regulation (Hancock et al., 2015). This approach enables us to reduce our two-state model (Equation (20)) to the following one-state model: assuming that disturbance is bounded and that 'weak' rapid buffering occurs, where b x is sufficiently large such that where Equation (22) can be related to buffering speed using We refer to Equation (22) as weak rapid buffering due to the extra condition for a bounded disturbance. This allows for a standard vector field definition and non-dimensionalization using a general estimate of the largest typical Dy.
We refer to 1/c = b x +g x as the buffering speed. Since b y = B/c, we can describe buffering parameters in terms of the buffering equilibrium ratio B and buffering speed 1/c, instead of b x and b y .
For the case in the main text, where g x = 0, then Equation (21)  where the feedback gain can be written In our analysis, we treat H and B as free parameters in order to study their particular effect on regulation. Note our treatment of H as a free parameter despite its dependence on B (via G). In any scenario where B is understood as varying independently of H, it is therefore assumed that either h y varies correspondingly, in order to keep H constant, or that g x and g y vary accordingly, in order to keep G constant. H and B are always independent for the case g x = 0, as presented in the main text.
We can also note that for linear removal, G = y=p y , and so the feedback gain can be written in the form H = À y p y vp y ðyÞ vy : To derive the reduced model in Equation (21), we use the total quasi steady-state approach (Borghans et al., 1996). We first transform the two-state model in Equation (A.20) to where Dy T = Dy + Dx is the total concentration. Treating Dy T as the slow variable and Dx as the fast variable, we assume that Dx is at quasi-steady state (D _ xz0), and so D _ y T = À À g y + h y Á Dy À g x Dx + d y ðtÞ Dx = BDy: We note that the buffering equilibrium ratio at quasi-steady state is B = Dx Dy for non-zero Dy. Substituting, we have D _ y T = À À g y + Bg x + h y Á Dy + d y ðtÞ = À À G + h y Á Dy + d y ðtÞ: Applying Dy T = Dy + Dx to the slow time scale, we have and thus obtain Equation (21). To derive the time-scale separation condition in Equation (22), we need to place the two-state model into standard singular-perturbation form (Khalil, 2002). Eliminating Dy from Equation (23), we have the two-state model We next non-dimensionalize the equations with the aim of ensuring that all terms are of a similar scale, with the exception of a small parameter 3. We initially choose an arbitrary time scale, arbitrary maximum of Dy T , and scale Dx in terms of Dy T . We then select a typical maximum value of Dy T and a typical time scale.
Non-dimensionalizing the variables, we have such that t is the non-dimensionalized time scale and where the disturbance is non-dimensionalized to normalize the vector field. We next selectỹ T , T s , and T f based on a typical maximum value of Dy and d y , which we can use to transform the equations into a standard singular perturbation form. We have dDy Tn dt = À l 1 Dy Tn + ðl 2 À l 3 ÞDx n + d n ðtÞ 3 dDx n dt = Dy Tn À Dx n ; ( Equation 24) where and where T s and T f are the slow and fast time scales, respectively. The assumption in Equation (22) follows from writing 3(1 in expanded form.

A.4. Steady-State Disturbance Rejection
The following section analyzes the effect of steady-state disturbances on the minimal model.
To understand the effect of disturbances, we first analyze steady-state disturbances. This is an important regulatory case, and also illustrates the different sensitivity functions used to characterize disturbance rejection.
To obtain the steady-state deviation from the nominal under a constant and persistent disturbance d y = d y , we set D _ y = D _ x = 0 in Equation (20). This yields where d y is a constant. It is important to note that G is dependent upon B only if g x s0.
In order to evaluate the effectiveness of regulation strategies for steady-state disturbance rejection, we introduce a sensitivity function, defined for two separate cases. One for the case where disturbances scale proportionally to p y , the steady-state production of y: and another for the case where disturbances are independent from other signals: In both expressions, jDyj is normalized by y, its nominal steady state. In the former case, we also normalize the disturbance, d y , by p y , the nominal production rate (with feedback) of y. This allows us to compare two systems that achieve the same nominal y using different rates of production p y . For example, the addition of a buffer to a system requires, as compensation, an increase in the rate of production of y or a decrease in its effective rate of removal in order to maintain the same y. In the case of an increased nominal production rate, disturbances that enter the system via y production will be amplified. Note that other scenarios, such as when disturbances scale with the rate of removal of y, are not analyzed here. Using Equations (19) and (26), we can simplify Equations (27) and (28)  For the case of independent disturbances, increasing the amount of buffer typically increases the total removal rate as both species are degraded, and so in steady state, the total production also increases to keep the steady state at the required concentration. This increase in production makes the relative effect of a perturbation smaller. However, this form of regulation can also be achieved by any increase in production/removal, without requiring the need of a buffer.

A.5. Oscillatory Disturbance Rejection: Rapid Buffering
The following section analyzes the effects of oscillatory disturbances on the minimal model, which generalizes the analysis in the previous section. Fourier (and the related Laplace) transforms are used for analysis. Fourier transforms allow a disturbance to be decomposed into its constituent oscillatory signals (frequencies), and for the effect of each constituent frequency on the system to be quantified.
In this section, we determine the sensitivity of the class of rapid buffering systems to oscillating disturbances of different frequencies. Oscillating disturbances with zero frequency are the special case of a constant disturbance analyzed in A.4. For dynamic disturbances, we define the sensitivity function to be where the square-root of the power of a signal is defined as Dy 2 dt s : The power and corresponding sensitivity function are useful for analyzing persistent disturbances, i.e., those that do not fade away over time.
To determine the sensitivity to oscillating disturbances, we apply a frequency domain approach with Laplace/Fourier transforms. From Equation (21), the model for rapid buffering is ð1 + BÞD _ y = À Gð1 + HÞDy + d y ðtÞ; which has a frequency domain representation of ð1 + BÞsYðsÞ = À Gð1 + HÞYðsÞ + D y ðsÞ; where Y(s) = L{Dy(t)} is the output and D y (s) = L{d y (t)} is the disturbance in the frequency domain, L(,) is the Laplace transform, and s is the complex frequency. The transfer function from the disturbance D y (s) to the output Y(s) is GðsÞ = YðsÞ D y ðsÞ = 1 ð1 + BÞs + Gð1 + HÞ : Setting s = iu (where i is the imaginary unit number), we have (see also Figure S2) where G(iu) denotes the frequency response from D y (iu) to Y(iu). The magnitude of this frequency response is given by GðiuÞ 2 = 1 ð1 + BÞ 2 u 2 + G 2 ð1 + HÞ 2 : (Equation 30) Using the properties of the transfer function (Parseval's Theorem for the equivalence of the power of a signal in the time or the frequency domain (Doyle et al., 1990)) and Equation (19), we obtain Normalizing the disturbance frequency u by the unregulated system time scale T u = 1/g y , we have where u n = uT u . If we assume that g x = 0 and that the removal is linear, then G = g y + Bg x = g y and G = G = g y . Thus, for the case presented in the main text we have f u = 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ð1 + BÞ 2 u 2 n + ð1 + HÞ 2 q :

(Equation 31)
A.6. Impulse Disturbance Rejection: Rapid Buffering The following section analyzes the effect of impulse disturbances on the minimal model. The analysis is completed using Fourier transforms, allowing the overall effect to be calculated by summing (via integration) the effect from all the constituent signals of the frequency-decomposed impulse disturbance. We next determine the sensitivity of the model in Equation (18) to an impulse with rapid buffering. The sensitivity function f imp used for analysis is where T u = 1/g y is the time scale of the unregulated system, b d y is a measure of the impulse disturbance, and the L 2 norm is defined as which is a metric of the accumulated deviation of Dy over time. The impulse disturbance is d y ðtÞ = b d y dðtÞ where d(t) is a unit impulse and b d y is a constant. The unit impulse is (informally) defined to be d(t) = N when t = 0 and d(t) = 0 when ts0, such that R N ÀN dðtÞdt = 1 (Skogestad and Postlethwaite, 2005). The non-dimensionalization of kDyk L2 by ffiffiffiffiffiffiffiffiffi ffi is based on y 2 corresponding to the integrand Dy 2 in Equation (33) and T u corresponding to the integration over time. The non-dimensionalization of d y (t) is by p y T u , where p y corresponds to the integrand and T u is due to the integration over time.
For linear removal with g x = 0, as presented in the main text, then G = g y + Bg x = g y and so G = g y . This gives p y T u = p y =g y = y, and the sensitivity function reduces to For an impulse d y = b d y dðtÞ, we cannot use the 'weak rapid buffering' assumption given by Equation (22), as the impulse does not have a finite maximum. Instead, we provide the alternate 'strong rapid buffering' assumption For the general case (g x s0), the sensitivity function for strong rapid buffering (Equation (34)) becomes where G is the steady-state removal rate, T u = 1/g y is the unregulated open-loop timescale, and G = g y + Bg x is the local removal rate.
For the case where g x = 0 and removal is linear, G = g y + Bg x = g y and G = G = g y . For this case, we also note that p y T u = p y =g y = y. Thus, for the case presented in the main text, we have A.7. Impulse Disturbance Rejection: Non-Rapid Buffering The following section carries out further analysis on the effect of impulse disturbances on the minimal model. The analysis is similar to Section A.6, but with a higher-order model. We next characterize impulse disturbance rejection when buffering is no longer assumed to be rapid. For this, we determine the frequency response associated with the linearized two-state minimal model The frequency domain representation is sYðsÞ = À b y YðsÞ + b x XðsÞ À g y YðsÞ À h y YðsÞ + D y ðsÞ sXðsÞ = b y YðsÞ À b x XðsÞ À g x XðsÞ; which leads to where the transfer function from D y (s) to Y(s) is e7 Cell Systems 5, 1-11.e1-e23, November 22, 2017 (Equation 36) where 1/c = b x +g x is the buffering speed.
To calculate the impulse-response sensitivity function for the general case, we need to determine the 2-norm of G(s  (36) and (37), we have From Equation (32), the sensitivity function is and T u = 1/g y . If k/1, then the non-rapid sensitivity function approaches the rapid case in Equation (35). For k z 1, we require that which is an alternative form of the strong buffering assumption given by Equation (34). The strong rapid buffering condition in Equation (34) requires a larger b x than the weak rapid buffering condition given in Equation (22). We cannot use Equation (22) as impulses do not have a maximum, and so do not fit into the standard singular-perturbation form in Equation (24). Writing Equation (39) in terms of normalized buffering speed R instead of c, we have For the non-rapid case where removal is linear and g x = 0, resulting in G = G = g y + Bg x = g y , we have

A.8. Optimal Impulse Disturbance Rejection
The following section carries out further analysis on the effect of impulse disturbances on the minimal model. Specifically, it uses the results of the previous section to minimize f imp and show that the optimal regulation occurs when the buffering speed is extremely rapid. We next show that strong rapid buffering (Equation (34)) is optimal for rejection of an impulse disturbance at d y . In Equation (39), changing the buffering speed only changes k. Following from Equation (34), the strong rapid buffering limit is therefore, for fixed B and H and removal constants, if we take the limit of strong rapid buffering, then k/1 and f imp in Equation (39) approaches its infimum value.

A.9. Stability Analysis with Feedback Delays
The following section analyzes the stability of the minimal models. To do so, the models of the previous sections are extended to incorporate a feedback delay. Stability conditions are then determined using the Nyquist criterion (Nise, 2004)-a method using Fourier/ Laplace transforms of the open-loop system (i.e., with disconnected feedback) to determine the properties of the closed-loop system (i.e., with connected feedback).
In this section, we analyze the stability of a class of systems with feedback delay and show that buffering can help stabilize an otherwise unstable feedback system. We carry out this analysis for both rapid and non-rapid buffering. Stability Analysis: Rapid Buffering. We first complete stability analysis under the assumption that buffering is rapid. With rapid buffering, the model with a feedback delay is a modified version of Equation (21) The closed-loop transfer function is thus given by (see also Figure S3) ( Equation 41) To analyze the class of system with delays, we use the Nyquist stability criterion and the open-loop transfer function (Nise, 2004). The open-loop transfer function from U f (s) to the output Y(s) is given by (including buffering but not feedback and considering d y (t) = 0 in Equation (40)) For the stability of this case, the Nyquist stability criterion requires that h y GðiuÞ <1 at the frequency u = u c where :Gðiu c Þ = p.
Starting with the characterization of :Gðiu c Þ = À p, we have :Gðiu c Þ = À u c t À tan À1 ðu c TÞ = À p; where T = 1 + B G . If we asymptotically approximate tan À1 ð,Þ using tan À1 ðu c TÞz we obtain :Gðiu c Þz 8 > < > : : This yields The Nyquist stability criterion (Nise, 2004)  which is equivalent to resulting in the constraint where f u,low represents f u for low frequencies, i.e., for values of u where (1 + B)u((1 + H) in Equation (31). Thus, the feedback system is stable if and only if the feedback gain, H, is smaller than a specific upper bound (given in Equation (44)) that depends on the amount of delay in the feedback loop (u c depends on t in Equation (42)). From Equation (42), we can write where Cell For the case of g x = 0 presented in the main text, Equations (43) and (44)  s ; a = max½p=2; pt=ðt + TÞ; and T = ð1 + BÞ g y : Graphical Stability Analysis: Non-Rapid Buffering. We next look at how stability conditions are changed for the case where buffering is not rapid. To do so, we derive the transfer function for the two-state model in Equation (20), modified to include feedback delay. We then graphically show an increase in stability gain margin as buffering is slowed from rapid buffering. Starting from the initial two-state model in Equation (20), the two-state model with delay can be written as D _ y = u f ðt À tÞ À b y Dy + b x Dx À g y Dy + d y ðtÞ D _ x = b y Dy À b x Dx À g x Dx u f ðtÞ = À h y DyðtÞ: Its Laplace Transform is given by sYðsÞ = e Àst U f ðsÞ À b y YðsÞ + b x XðsÞ À g y YðsÞ + nðsÞ sXðsÞ = b y YðsÞ À b x XðsÞ À g x XðsÞ U f ðsÞ = À h y YðsÞ; and the open-loop transfer function from U f (s) to Y(s) is where 1/c = b x +g x is the buffering speed. In Figure S4, we can see an example of the increase in gain margin from slowing buffering from the rapid limit. The increase in gain margin occurs due to the increase in phase of the open-loop system.

A.10. Molecular Noise Analysis
The following section analyzes the molecular noise in the minimal models. The system of chemical reactions is simulated for the case where the variables represent discrete numbers of molecules. We use large volume/numbers of molecules and linearization approximations to derive a linear Chemical Langevin Equation, in which variables are continuous concentrations, from the biochemical reactions. This simplified equation is then analyzed to determine the magnitude and frequency composition of the noise at the regulated species.
In this section, we determine the molecular noise due to the discrete molecular nature of the reactions. For simplicity and to better complement preceding sections, we use Chemical Langevin Equations (CLE) with the linear noise approximation to analyze molecular noise with linear removal terms. The CLE with the linear noise approximation is a larger-volume, linearization approximation of the Chemical Master Equation (Klipp et al., 2009;Bressloff, 2014;Elf and Ehrenberg, 2003)-equations representing the stochastic system for discrete molecular levels. The linear noise approximation also exactly agrees with the Chemical Master Equation up to second-order moments for chemical systems composed of zero, first, and some second order reactions (Grima, 2015). We use the Gillespie algorithm (Klipp et al., 2009) for simulations of discrete molecular levels presented in the main text.
If we assume linear removal of species y and ignore removal of species x, the biochemical system has the transitions y n ! p yn ðy n Þ y n + 1 y n ! g y y n y n À 1 ðy n ; x n Þ! g yn ðy n ; x n Þ ðy n À 1; x n + 1Þ ðy n ; x n Þ! g xn ðx n ; y n Þ ðy n + 1; x n À 1Þ; where y n = Uy is the number of molecules of y in a system with volume U, x n = Ux is the number of molecules of x, p yn (y n ) = Up y (y n /U), g yn (y n ,x n ) = Ug y (y n /U,x n /U), and g xn (x n ,y n ) = Ug x (x n /U,y n /U). Linear buffering rate equations are used in the simulations. We next use a linear noise approximation to determine the stationary variance. The CLE used for analysis with large volume/ numbers of molecules is (Klipp et al., 2009;Bressloff, 2014) dy dt = p y ðyÞ À g y ðy; xÞ + g x ðx; yÞ + b x x À g y y + U À1=2 ffiffiffiffiffiffiffiffiffiffiffi p y ðyÞ q x 1 À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi g y ðy; xÞ q x 2 + ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi g x ðy; xÞ p x 3 À ffiffiffiffiffiffiffi ffi g y y p x 4 dx dt = g y ðy; xÞ À g x ðx; yÞ + U À1=2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi g y ðy; xÞ where x i s are Gaussian white noise sources with zero mean and unit variance, p y includes feedback, and the equations are presented using a standard derivative notation for Langevin equations. We next use the linear noise approximation about the steady state of the deterministic system (Klipp et al., 2009;Bressloff, 2014;Elf and Ehrenberg, 2003), noting that the mean of each concentration is assumed to be equal to its deterministic steady state. We set y = y + U À1=2 b y and where g y = g y ðy; xÞ, g x = g x ðx; yÞ, and the remaining coefficients are defined in Equation (20). When buffering is linear, then g y = b y y and g x = b x x. We can write Equation (46) where b z = ðb y; b xÞ T is the vector of states, x is the vector of Gaussian white noise sources, The variance of b z can be determined from the covariance matrix Q = covðb zÞ by solving (Klipp et al., 2009;Bressloff, 2014;Elf and Ehrenberg, 2003) b and the covariance matrix Q n of (y n , x n ) can be determined by solving (Oyarzú n et al., 2015;Elf and Ehrenberg, 2003) b Solving for the variance s 2 yn of y n , we have where g yn = Ug y , hy n i = Uhyi = Uy and The first term in Equation (47) represents noise from production and removal reactions, while the second term represents noise due to the buffering reactions. The resulting noise intensity is The function f noise is also known as the noise strength or coefficient of variation. The noise intensity for rapid buffering (r/0) and linear buffering (g y = b y y), which is presented and discussed in the main text, is It can be seen that for slower buffering (where r > 0), feedback can regulate the buffering molecular noise term. In this case, the production/removal noise can be attenuated by both increased feedback gain, H, and increased buffering equilibrium ratio, B. Also, non-linear buffering can reduce the contribution of buffering noise if g y <b y y, under linear noise approximation assumptions.
We next analyze the linearized equation using a frequency-domain approach, noting that the Laplace transform is used interchangeably with the Fourier transform by setting s = iu. The equivalence of a frequency-domain methodology can be seen by noting that (Doyle et al., 1990) Varðb Using the independence of white noise sources, the frequency-domain representation of Equation (46) where b Y is the frequency-domain representation of b y, N 1 to N 4 are the frequency-domain representations of x 1 to x 4 , respectively, and To analyze the noise bandwidth for rapid buffering, we set s = iu in Equation (49) We can observe that G y acts as a low-pass filter with a bandwidth between 0 and g y (1 + H)/ (1 + B), while G b acts as a band-pass filter with a bandwidth between g y (1 + H)/(1 + B) and (1 + B)/c. By noting that white noise has a uniform spectral density, we can see that the buffering reaction introduces higher-frequency noise while the production and removal reactions introduce lower-frequency noise. A.11. Comparison to Technological Feedback Control The following section compares buffer-feedback systems to technological controllers. Similarities and differences are shown using Laplace transforms, as this transform is commonly used for the design and representation of technological controllers.
In this section, we compare buffering to derivative feedback in the linear case. We first show that rapid buffering is mathematically equivalent to negative derivative feedback. This finding generalizes a previous in silico observation that creatine acts as derivative control for ATP (Cloutier and Wellstead, 2010). We also show the equivalence of (general) buffering to either derivative filtering or lead controllers, both closely related to pure derivative controllers. Finally, we represent the class of buffering-feedback systems in block diagram form and show the similarity between measurement noise and the molecular noise produced by buffering.
PID (Proportional-Integral-Derivative) feedback controllers provide a straightforward method for considering and combining past (integral), present (proportional), and future (derivative) state information (Aströ m and Murray, 2008). In particular, derivative feedback works by opposing and cancelling a portion of the would-be (i.e., short-term future) rate of change of the regulated species, and in a manner that is not dependent on current states (i.e., current concentrations). In our biochemical model, rates of change represent net fluxes. When d y and therefore the inflow rate of y increases, a buffer reduces the would-be rate of increase to the net flux of y by conversion to the buffer species x. Conversely, when d y decreases, a buffer reduces the would-be rate of decrease to the net flux of y by x/y conversion. The rapid equilibrium of the buffering reaction ensures that the rate of conversion between y and x depends only on the net flux of y and not its current concentration; this is because the concentration-dependent contributions to the y4x conversion rates, b y y and b x x, are always balanced.
To show the equivalence with derivative feedback, we use the model for rapid buffering in Equation (21) with g x = 0: Treating the buffer as a second 'feedback' term, this model can be rewritten as D _ y = À g y Dy À h y Dy À BD _ y + d y ; in which we observe that the buffer term, ÀBD _ y, is equivalent to a negative derivative feedback term with gain B. The model in Equation (50) is thus equivalent to a negative feedback model comprising a proportional feedback term Àh y Dy with proportional feedback gain h y , and a derivative feedback term ÀBD _ y with derivative feedback gain B. If we cannot assume rapid buffering or that g x = 0, we use the two-state model in Equation (20), which in the frequency domain can be written as sYðsÞ = À h y YðsÞ À b y YðsÞ + b x XðsÞ À g y YðsÞ + D y ðsÞ sXðsÞ = b y YðsÞ À b x XðsÞ À g x XðsÞ; where Y(s) and X(s) are the Laplace transforms of Dy and Dx, respectively. The second equation can be rearranged as and substituting this into the first equation gives us sYðsÞ = À h y YðsÞ À g y YðsÞ + b x b y s + b x + g x À b y YðsÞ + D y ðsÞ: Isolating the contribution of buffering into a feedback term called U b (s), we can rewrite the model as sYðsÞ = À h y YðsÞ À g y YðsÞ + U b ðsÞ + D y ðsÞ where C d (s) represents the buffering subsystem. This subsystem can be rewritten as which is the standard form of a lead controller (Nise, 2004). Thus for g x s0, buffering is equivalent to a lead controller, a form of controller closely related to derivative controllers (Nise, 2004).
Cell Systems 5, 1-11.e1-e23, November 22, 2017 e14 For the case where the disturbance is on the buffering species and feedback is on the regulated species ( Figure 1A in the main section of the paper with disturbance d x and not d y ), we consider the minimal model with nominal steady state ðy; xÞ and deviations Dy = y À y and Dx = x À x. Linearization of this model yields where h y = À v vy p y ðyÞj y = y ; b y = v vy À g y ðy; xÞ À g x ðx; yÞ Á j ðy; xÞ = ðy;xÞ ; b x = v vx À g x ðx; yÞ À g y ðy; xÞ Á j ðy;xÞ = ðy;xÞ ; Taking the Laplace transform gives us which can be rearranged as or furthermore as using Y(s)/D y (s) from Equation (36). Setting s = iu gives us Normalizing the disturbance frequency u by the unregulated system time scale T u = 1/g y , such that u n = uT u , we have For an oscillating disturbance, we have the sensitivity function For the case presented in the main text where g x = 0 and so q = 1, the sensitivity function becomes We next look at the impulse response. Calculating the norm of G x (s) using Equations (37) and (38), we have impulse. This leads to the sensitivity function From Equation (54), it can be seen that reducing R for a given B improves the rejection of impulse disturbances at d x .

A.13. Slow Buffering with Reservoir Feedback
The following section carries out analysis on the effect of feedback from the buffering species instead of the regulated species. The section shows that the feedback can act as a form of integral feedback for slow buffering. The analysis uses Laplace tranforms to show the similarities to proportional and integral control, in a similar approach to Section A.11.
We consider the revised minimal model in which the feedback has been altered from p y (y) to p y (x). Assuming equivalent steady states to Equation (18), this model can be linearized as e17 Cell Systems 5, 1-11.e1-e23, November 22, 2017 where h x = v vy p y ðxÞj x = x . Written in the frequency domain, we have sYðsÞ = À h x XðsÞ À b y YðsÞ + b x XðsÞ À g y YðsÞ + D y ðsÞ sXðsÞ = b y YðsÞ À b x XðsÞ À g x XðsÞ; which leads to Thus, the feedback term can be written as YðsÞ: For rapid buffering, we have the feedback approximation which is the similar to feedback from y, only differing by the gain. We observe that there is a pole at s p = À(b x +g x ), where the transfer function is infinite. For the slow buffering limit, the pole s p = À(b x +g x )/0, and so we have the feedback approximation which is a form of integral feedback.
Combining the effect of both buffering and feedback from the buffering species, we have which can take the more general lead or lag controller form, rather than just being limited to a lead controller. For the slow buffering limit, where the pole s p = À ðb x + g x Þ/0, we have an approximate PI (proportional plus integral) regulator

B. Regulation of Glycolysis
This supplementary section introduces a minimal model for glycolysis, analyzes the model, and uses the qualitative results to find new relationships in existing experimental data.
In this section we examine the regulation of glycolysis. In particular, we look at the effects of glycolytic feedback and creatine phosphate buffering on the stability and regulation of ATP concentrations. We then use this analysis to compare the regulation of ATP by different species.
The minimal models of glycolysis that we analyze here are built on models found in existing literature that do not include buffering (Chandra et al., 2011). Initially, we look at the case where buffering is rapid (B.1), then analyze more general cases of buffering (B.2), followed by cases where more general feedback structures are allowed (B.3). As in (Chandra et al., 2011), we assume that the total concentration of adenosine phosphates in the cell (ATP + adenosine diphosphate (ADP) + adenosine monophosphate (AMP)) remains constant, and that the activating effects of AMP can be modelled as ATP inhibition. By further assuming total creatine (Cr + pCr) to be constant, we simplify the buffering reaction as pCr4ATP. where y represents ATP concentration, x represents pCr concentration (for vertebrates), and m is a lumped variable representing intermediate metabolites (see also Figure 5F). Also, g y (y,x) and g x (x,y) are the reaction rates for the ATP 4 pCr reaction, the constant V y is the nominal ATP demand, and q is the stoichiometric investment/return ratio for the autocatalysis of ATP.
The reaction rates can, for example, be modelled by Hill functions, such that fðyÞ = V PFK y a 1 + g h y h ; wðyÞ = V PK 1 + g g y g ; g y ðy; xÞ = b y y where V PFK , g h , h, and a are parameters describing the reaction rate of phosphofructo kinase (PFK), and V PK , g g , g are the parameters describing the reaction rate of pyruvate kinase (PK). The parameters b y , b x , K y and K x describe the reaction rates of creatine kinase, and are also dependent upon creatine and ADP concentrations. The steady-state solution of Equation (55) is where B is the steady-state ratio between ATP and pCr, such that g y ðBx; xÞ = g x ðx; BxÞ.

B.1. Rapid Buffering
The following section analyzes a (further simplified) minimal model for glycolysis, determining the critical feedback gain for which its steady state is stable. The stability analysis is carried out by linearizing the system about its steady state, followed by eigenvalue analysis of the linear system. The critical feedback gain is then used to determine a limit on the effectiveness of steady-state regulation.
We first look at the case where buffering is rapid. Linearizing Equation (55) about the nominal steady state with deviations Dy = y À y, Dm = m À m, and Dx = x À x, results in a f = vf vy ðyÞ; a w = vw vy ðyÞ: The buffering equilibrium ratio is For the rapid case, we can determine a reduced model similar to that used in STAR Methods A.3. By assuming rapid equilibrium, Dx = BDy, and that the total concentration Dy T = Dy+Dx is a slow variable, the model in Equation (57)  where ÀwðyÞ ð a f À a w mÞ y V y For stability we require the trace Tr(A 1 )<0 and determinant Det(A 1 )>0, such that ÀwðyÞ À 1 1 + B ðqa f À ðq + 1Þa w mÞ < 0 and ðqa f À ðq + 1Þa w mÞ 1 1 + B wðyÞ À ða f À a w mÞðq + 1Þ 1 1 + B wðyÞ > 0; which simplifies to Àqa f < ð1 + BÞwðyÞ À ðq + 1Þa w m and À a f > 0: In non-dimensionalized form, this can be rewritten as a constraint on feedback gain: The steady-state regulation is lower bounded by When the reaction rates are modelled by the Hill functions in Equation (56), we have H PFK = h g h y h 1 + g h y h À a and H PK = g g g y g 1 + g g y g :

B.2. Non-Rapid Buffering
The following section carries out further analysis of a minimal model of glycolysis. The analysis is similar to Section B.1, but with a higher-order model to represent non-rapid buffering.
We next study stability constraints for the general buffering case. For simplicity, we first determine stability conditions before nondimensionalizing. For the non-rapid case, Equation (57) can be rewritten as 2 has the characteristic polynomial s 3 À ða + d + gÞs 2 + ðad + dg + ag À ef À bcÞs + ð À adg + aef + bcgÞ = 0: Cell Systems 5, 1-11.e1-e23, November 22, 2017 e20 Stability requires that ah À ða + d + gÞ > 0 bhad + dg + ag À ef À bc > 0 gh À adg + aef + bcg > 0 ab > g: Therefore, for the buffering model, stability requires that and, from ab>g, that When non-dimensionalized, the stability conditions can be written as where H PFK = À a f y V y = À vf vy ðyÞ y fðyÞ H PK = À a w m y V y = À vw vy ðyÞ y wðyÞ We can check special cases of the general stability conditions by taking limit cases of the buffering speed b x . As b x /N, the stability conditions are which is identical to the 'rapid' buffering case. As b x /0, the stability conditions are which is equivalent to the case without buffering.
The classical loop transfer function is composed of two components corresponding the feedback and buffering, such that L = L B + L H , where L B = Bsðs + wÞ s 2 + ½w + qa f0 À ðq + 1Þa w ms À wa f0 L H = C h ðw À qsÞ s 2 + ½w + qa f0 À ðq + 1Þa w ms À wa f0 : In (Chandra et al., 2011), it was shown that the right-half plane (RHP) poles and zeros of the classical loop function fundamentally limits the regulation of glycolysis. It can be seen that PFK feedback is associated with a RHP zero at z = w/q. However, buffering is associated with a left-half plane (LHP) pole at z = Àw, and so is not subject to the same fundamental constraints.
For proportional feedback, L has zeros at If this case holds, there are no fundamental constraints on feedback due to RHP zeros. Biologically, the RHP zero corresponds to the initial drop in ATP due to the investment stage of glycolysis, during which the ATP concentration initially decreases. The shifting of the zero from the RHP to the LHP corresponds to the buffer supplying the initial ATP investment to PFK.
We next assume that feedback may take any stable, causal form. L(s) has zeros at C h ðsÞðw À qsÞ + Bsðs + wÞ = 0: We can write C h ðsÞ = h ðs À z 1 Þ.ðs À z m Þ ðs À p 1 Þ.ðs À p n Þ ; where m % n, and so sðs + wÞðs À p 1 Þ.ðs À p n Þ À qh B ðs À z 1 Þ.ðs À z m Þ s À w q = 0: (Equation 61) Equation (61) is a 'reverse' root locus problem if qh/B is treated as the gain, noting that here we are using the method to study the open-loop zeros rather than the closed-loop poles. The zeros of L(s) start at 0, w, and the poles of the feedback. As qh/B is increased from zero, the zeros of L(s) are all in the LHP, noting that the zero starting at z = 0 moves to the left (left of an even number of poles/ zeros on the axis for 'reverse' root locus). Thus, for small qh/B, L(s) has no RHP zeros. Or, in other words, the fundamental limits associated with RHP zeros are removed for a sufficiently large buffering equilibrium B. B.4. Discussion: Comparison to Experimental Data The following section compiles published experimental data for glycolytic feedback vs buffering within several different organisms. The published data shows a positive trend in feedback gain vs buffering equilibrium ratio. This trend, along with observed cases of oscillations that occur when buffering is removed, supports the notion of a critical feedback gain that is dependent upon buffering, as presented in the last section.
Motivated by the mathematical analysis, we investigated for an empirical relationship between glycolytic feedback via PFK and ATP buffering by pCr. Specifically, we compared proxies for the feedback gain, H PFK , and the buffering equilibrium ratio, B, of different species. The feedback gain from ATP concentration to PFK activity served as a proxy for H PFK , while the steady-state ratio B = [pCr]/[ATP] served as a proxy for the buffering equilibrium ratio B.
The feedback gain was estimated using H PFK = À y fðyÞ Df Dy ; where y is the nominal ATP concentration, fðyÞ is the nominal PFK activity, and Df/Dy is the slope of PFK activity to ATP concentration at the nominal point. Experimental values for these three quantities were determined from literature, directly from plots of ATP concentration vs PFK activity. Nominal values of ATP concentration and PFK activity were chosen at half maximum PFK activity, while the slope was estimated graphically at this same point. Lines of best fit were used to estimate the slope when available. The point of half maximum PFK activity was selected to estimate the maximum slope, which is the point most likely to cause instability. Using the point of half maximum PFK activity allowed a consistent approach across all data sets.
The data values that were selected from literature are presented in the table in Figure S7 and plotted in Figure 6. For further support of the results, we also used data from published studies that have observed glycolytic oscillations after the removal of creatine phosphate: in vitro, in cell-free extracts of rat skeletal muscle (Tornheim and Lowenstein, 1973;Tornheim et al., 1991), and in vivo, in rabbit ventricular myocytes (when aerobic metabolism is also inhibited) (Yang et al., 2008). These cases show that removal of ATP buffering can lead to feedback instability.
The stability boundary in Figure 6 is a qualitative approximation that does not take into account q or H PK in Equation (60). It is included for illustrative purposes to highlight that the removal of buffering can lead to instability. In Figure 6, we have, in the absence of experimental data, approximated the feedback gain in rabbit ventricular myocytes to be equal to that in rabbit skeletal myocytes, since both have predominantly the same PFK isoform. ν y (y) = γ y (y −ȳ) + ν y (ȳ) Figure S1: Example of a removal rate function in Equation (18), illustrating the significance of the linearization constant γ y and the steady-state constantΓ (assuming ν x (x) = 0). γ y is the slope of the removal rate function linearized about the equilibrium point.Γ is the slope of the line connecting the equilibrium point to the origin. Frequency (rad/s) Figure S3: A Bode plot of the closed-loop transfer function given by Equation (41) with nondimensionalized time. Feedback gain H = 0 or 1.5, buffering equilibrium ratio B = 0 or 10, and delay τ = 2. The feedback plot shows a stable case near instability, where the magnitude response contains a 'resonance spike'. The plot also shows that buffering removes the resonance spike and highlights the requirement for buffering to reject high-frequency disturbances rather than feedback.  Figure S6: The buffering-feedback regulation system when γ x = 0, presented in block-diagram representation. Y is the frequency transform of set point deviation ∆y.