††thanks<EMAIL_ADDRESS>
# Thermodynamic optimality of glycolytic oscillations

Pureun Kim  Korea Institute for Advanced Study, Seoul 02455, Korea  Changbong
Hyeon  Korea Institute for Advanced Study, Seoul 02455, Korea

###### Abstract

Temporal order in living matters reflects the self-organizing nature of
dynamical processes driven out of thermodynamic equilibrium. Because of
functional reason, the period of a biochemical oscillation must be tuned to a
specific value with precision; however, according to the thermodynamic
uncertainty relation (TUR), the precision of oscillatory period is constrained
by the thermodynamic cost of generating it. After reviewing the basics of
chemical oscillations using the Brusselator as a model system, we study the
glycolytic oscillation generated by octameric phosphofructokinase (PFK), which
is known to display a period of several minutes. By exploring the phase space
of glycolytic oscillations, we find that the glycolytic oscillation under the
cellular condition is realized in a cost effective manner. Specifically, over
the biologically relevant range of parameter values of glycolysis and
octameric PFK, the entropy production from the glycolytic oscillation is
minimal when the oscillation period is (5 – 10) minutes. Further, the
glycolytic oscillation is found at work near the phase boundary of limit
cycles, suggesting that a moderate increase of glucose injection rate leads to
the loss of oscillatory dynamics, which is reminiscent of the loss of
pulsatile insulin release resulting from elevated blood glucose level.

From the perspective of thermodynamics on the macroscale, it is clear that
living organisms are sustained far from thermodynamic equilibrium in that
energy and material currents are supplied to the systems and dissipated to the
environment Bustamante  _et al._ (2005). For human consisting of $\sim 37$
trillions of cells, the energy consumption is $\sim 100$ W. At the cellular
scale, this amounts to as much as $\sim 10^{9}$ $k_{B}T$ of energy being
consumed per second, which drives a vast number of cellular processes. Among
them, biochemical oscillations, exemplified with circadian rhythm Rust  _et
al._ (2011, 2007) and the cell cycle Novak and Tyson (1993); Tyson and Novak
(2001); Ferrell  _et al._ (2011), are the temporal dissipative structures that
emerge from molecular components and their dynamics in nonequilibrium
Prigogine (1978); Goldbeter (2018).

The full catabolism of glucose under aerobic condition is the primary source
of ATP; one glucose molecule could yield as many as 30-32 ATP molecules.
Besides its utmost importance along the metabolic pathway, the sustained
oscillation of substrate and product concentrations with time period of $\sim$
(5 – 10) min observed in the glycolysis of yeast, muscle, and pancreatic
$\beta$ cell is a phenomenon of great interest on its own Duysens and Amesz
(1957); Boiteux  _et al._ (1975); Goldbeter and Caplan (1976); Tornheim
(1997); Bertram  _et al._ (2004). Studies carried out on yeast indicate that
among a host of glycolytic enzymes, the allostery of phosphofructokinase
(PFK), which uses ATP hydrolysis to catalyze the phosphorylation of
fructose-6-phosphate (substrate, $S$) into fructose-1,6-biopohosphate
(product, $P$), was identified to lie at the core of the endogenous
oscillatory mechanism Hess and Boiteux (1971); Boiteux  _et al._ (1975). The
binding of product to the allosteric site induces a cooperative transition of
the PFK enzyme from its inactive to active state, introducing a nonlinear
positive response to the biochemical circuit. The product-activated allosteric
regulation of PFK enzyme becomes the source of the glycolytic oscillation.

In this paper, we study the phosphofructokinase (PFK) model for glycolytic
oscillations, and discuss how optimal this energy consuming process is in
achieving its desired oscillation period at the expense of thermodynamic cost,
which we quantify using the thermodynamic uncertainty relation (TUR), a
recently discovered thermodynamic principle Barato and Seifert (2015);
Gingrich  _et al._ (2016); Horowitz and Gingrich (2019); Hasegawa and Van Vu
(2019). In the section _Methods_ , we briefly review the TUR for dynamical
processes displaying temporal oscillations. Before studying the PFK model of
glycolytic oscillation, we review the Brusselator, a model chemical processes
exhibiting oscillations, and analyze it in detail at varying parameter values
by means of the phase diagram to identify the optimal condition of the process
in light of TUR. The _Results_ section covers the allosteric model of
octameric PFK, dynamic phase diagram, and the mass balance equations of PFK-
catalyzed substrate and product giving rise to the glycolytic oscillations. We
assess the optimality of glycolytic oscillations realized by octameric PFK in
terms of the period and fluctuations of the oscillations and the associated
thermodynamic cost of generating such dynamics, which can be quantified using
TUR. Finally, we will underscore the optimality of glycolytic oscillations
emerging from octameric PFK by comparing it with the dynamics anticipated for
hypothetical non-octameric PFK model.

##  I Methods

###  I.1 Thermodynamic uncertainty relation for oscillatory processes

Biological processes are inherently stochastic due to fluctuations in the
cellular environment. Suppose that time trajectories of a dynamical process
generated from a constant thermodynamic drive are available. Then, in order to
gain knowledge about the dynamics from a given trajectory with higher
precision, it is required to analyze a longer trajectory; however, generation
of a longer time trajectory incurs larger thermodynamic cost. According to the
thermodynamic uncertainty relation (TUR), there are trade-offs between the
total entropy production $\Delta S_{\text{tot}}(t)$ from the process generated
under a constant drive and the precision of a time integrated _current-like_
observable $X(t)$ to probe the dynamical process generated in nonequilibrium.
Furthermore, the product between $\Delta S_{\text{tot}}(t)$ and the squared
relative uncertainty, $\text{Var}[X(t)]/\langle X(t)\rangle^{2}$, defined as
the _uncertainty product_ $\mathcal{Q}$, is bounded below by twice the
Boltzmann constant ($2k_{\text{B}}$):

$\mathcal{Q}\equiv\Delta
S_{\text{\text{tot}}}(t)\frac{\text{Var}[X(t)]}{\langle X(t)\rangle^{2}}\geq
2k_{\text{B}}.$ (1)

Because of the lower bound, there is a minimal thermodynamic (entropic) cost
to generate a process with a certain precision, or that the precision of the
process is bounded by the thermodynamic cost being expended. The relation
applies to current-like output observables, with odd parity under time
reversal satisfying $X(t)=-X(-t)$ Hasegawa and Van Vu (2019), that arise from
dynamic processes that can be modeled using continuous time Markovian dynamics
in discrete network or in continuous space under _time-independent_ driving.
There have been further generalizations of the relation Proesmans and Van den
Broeck (2017); Horowitz and Gingrich (2019); Hasegawa and Van Vu (2019); Lee
_et al._ (2019); Agarwalla and Segal (2018); Lee  _et al._ (2018); Potts and
Samuelsson (2019); Koyuk  _et al._ (2018) including the those under time-
dependent driving Koyuk  _et al._ (2018); Koyuk and Seifert (2019, 2020).
Here, we confine ourselves to the original relation (Eq.1) since the
oscillatory dynamics considered in this study are generated under time-
independent driving. Throughout the paper, we set the Boltzmann constant
$k_{B}$ to unity for simplicity of the expressions.

There have been several applications of TUR to characterize biological
processes Pietzonka  _et al._ (2016); Hwang and Hyeon (2018); Uhl and Seifert
(2018); Song and Hyeon (2020); Piñeros and Tlusty (2020); Marsland III  _et
al._ (2019); Song and Hyeon (2021). For biological motors that move along
filaments, the net displacement of a motor $\Delta x(t)=x(t)-x(0)=-\Delta
x(-t)=x(0)-x(-t)$ satisfies the odd parity, with time translation symmetry,
and may be selected as a proper output observable to probe the functional
dynamics of the motor. For dynamical processes displaying temporal
oscillations, the time duration of the process, $t$, can be employed as a
current-like observable to monitor the progress of $n$ biochemical
oscillations, i.e., $X(t)=t$. In fact, it has recently been shown through the
large-deviation theory that in the limit of long times and large currents TUR
holds for the first passage times Gingrich and Horowitz (2017) as well, so
that a dynamical process consisting of $n(\gg 1)$ oscillations that has
occurred for a sufficiently long time duration $t$ satisfies the following
relation:

$\displaystyle\mathcal{Q}=\Delta S_{\rm tot}(t)\frac{\langle\delta
t^{2}\rangle}{\langle t\rangle^{2}}\geq 2k_{B}.$ (2)

The time duration, $t$, is decomposed as $t=\sum_{i=1}^{n}T_{i}$ where $n$ is
the number of oscillations for time $t$ and $T_{i}$ is a stochastic variable
representing the period of $i$-th oscillation. Then, the first and second
moments of the time duration are written as

$\displaystyle\langle t\rangle=\sum_{i=1}^{n}\langle T_{i}\rangle=n\langle
T\rangle$ (3)

and

$\displaystyle\langle t^{2}\rangle$ $\displaystyle=\sum_{i=1}^{n}\langle
T_{i}^{2}\rangle+2\sum_{i>j}\langle T_{i}T_{j}\rangle$ $\displaystyle=n\langle
T^{2}\rangle+n(n-1)\langle T_{i}T_{j}\rangle$ $\displaystyle=n\langle
T^{2}\rangle-n\langle T\rangle^{2}+\langle t\rangle^{2}$ (4)

where $\langle T\rangle\equiv(1/n)\sum_{i=1}^{n}T_{i}$ is the mean oscillatory
period. For renewal processes, the independence of two oscillatory periods,
$\langle T_{i}T_{j}\rangle=\langle T_{i}\rangle\langle T_{j}\rangle=\langle
T\rangle^{2}$, can be used, giving rise to the last line of Eq.4. Marsland III
_et al._ (2019) From Eqs. 3 and 4, one obtains $\langle t^{2}\rangle-\langle
t\rangle^{2}=n\langle T^{2}\rangle-n\langle T\rangle^{2}=n\times{\rm Var}(T)$,
which indeed holds for the Brusselator studied here (see Fig. S1). Even in the
presence of a finite correlation between oscillatory dynamics, $\langle\delta
T_{i}\delta T_{j}\rangle\simeq 0$ holds as long as the time gap corresponding
to $|i-j|$ is greater than the correlation time. Alternatively, in 1990s
Schnitzer and Block Schnitzer and Block (1995) showed for processive enzymatic
processes that the number of net catalytic events that have occurred for time
$t$, $n(t)$, is related to the enzymatic cycle time $\tau$, which corresponds
to the oscillatory period $T$ in this work, namely $\tau=T$, as
$\lim_{t\rightarrow\infty}\langle\delta n(t)^{2}\rangle/\langle
n(t)\rangle=\langle\delta\tau^{2}\rangle/\langle\tau\rangle^{2}=\langle\delta
T^{2}\rangle/\langle T\rangle^{2}$.

Taken together, the uncertainty product $\mathcal{Q}$ for processes
demonstrating temporal oscillations as follows Marsland III  _et al._ (2019);
Cao  _et al._ (2015); Morelli and Jülicher (2007):

$\displaystyle\mathcal{Q}$ $\displaystyle=\Delta
S_{\text{tot}}(t)\frac{\langle t^{2}\rangle-\langle t\rangle^{2}}{\langle
t\rangle^{2}}$ $\displaystyle\simeq\left[\frac{\Delta
S_{\text{tot}}(t)}{n\langle T\rangle}\right]\frac{\text{Var}(T)}{\langle
T\rangle}$ $\displaystyle=\dot{S}_{\text{tot}}\frac{\text{Var}(T)}{\langle
T\rangle}$ $\displaystyle=\Delta S_{\text{cyc}}\frac{\text{Var}(T)}{\langle
T\rangle^{2}}\geq 2k_{B},$ (5)

where $\text{Var}(T)=\langle(\delta T)^{2}\rangle=\langle T^{2}\rangle-\langle
T\rangle^{2}$. $\dot{S}_{\text{tot}}=\Delta S_{\text{tot}}(t)/n\langle
T\rangle=\Delta S_{\text{tot}}(t)/t$ is the rate of entropy production, and
the entropy production per cycle can be defined as $\Delta
S_{\text{cyc}}\equiv\Delta
S_{\text{tot}}(t)/n=\dot{S}_{\text{tot}}\times\langle T\rangle$. The
inequality in the last line of Eq.5 is identical to the expression
$\mathcal{N}\equiv\langle T\rangle^{2}/D\leq\Delta S_{\text{cyc}}/2$
($k_{B}=1$) of Marsland _et al._ ’s Marsland III  _et al._ (2019); Marsland
III and England (2017) where $D\equiv\text{Var}(T)=\langle(\delta
T)^{2}\rangle$ and $\mathcal{N}$ corresponds to the number _coherent_
oscillations. If the uncertainty product $\mathcal{Q}$ evaluated for a process
displaying temporal oscillations is close to the theoretical minimum $2k_{B}$,
one could presume that the process is close to its optimal condition where the
thermodynamic cost of generating an oscillatory dynamics with a certain
temporal precision is minimal.

For a given oscillatory process arising from a set of kinetic equations, it is
straightforward to evaluate the mean ($\langle T\rangle$) and variance of
oscillation period ($\text{Var}(T)$) in Eq.5 from the time trajectories of
dynamical variables, ${\bf x}(t)$, that satisfy ${\bf x}(t)\simeq{\bf
x}(t+T)$. To obtain the total entropy production rate from ${\bf x}(t)$, one
can utilize the evolution equation of probability density, $P({\bf x},t)$,
namely the Fokker-Planck (FP) equation. The FP equation is obtained from a
corresponding set of chemical master equations (CME) at a finite volume
($\Omega$) (see SI):

$\displaystyle{\frac{\partial P({\bf x},t)}{\partial t}}$
$\displaystyle=-\vec{\nabla}\cdot\left[{\bf H}\cdot\vec{F}({\bf x})-{\bf
D}({\bf x})\cdot\vec{\nabla}\right]P({\bf x},t)$
$\displaystyle=-\vec{\nabla}\cdot\vec{J}({\bf x},t),$ (6)

where ${\bf D}({\bf x})(=k_{B}T{\bf H}({\bf x}))$ is the $\Omega$-dependent
diffusion tensor with ${\bf H}({\bf x})$ being the motility tensor, and

$\displaystyle\vec{J}({\bf x},t)={\bf H}({\bf x})\cdot\vec{F}({\bf x})P({\bf
x},t)-{\bf D}({\bf x})\cdot\vec{\nabla}P({\bf x},t)$ (7)

is the probability current, $\vec{F}({\bf x},t)$ is the driving force vector.
The total entropy production rates, contributed by both the system and
reservoir, are obtained by averaging the corresponding _trajectory-based_
entropy production rates over the probability density $P({\bf x},t)$,

$\displaystyle\dot{S}_{\text{tot}}(t)=\langle\dot{s}_{\text{tot}}(t)\rangle=\int\dot{s}_{\text{tot}}(t)P({\bf
x}(t),t)d{\bf x},$ (8)

where
$\dot{s}_{\text{tot}}(t)=\dot{s}_{\text{sys}}(t)+\dot{s}_{\text{res}}(t)$ with
$s_{\text{sys}}(t)=-\log{P({\bf x}(t),t)}$ and
$\dot{s}_{\text{res}}(t)(=\vec{F}({\bf x},t)\cdot\dot{\bf x})/T_{r}$ where
$T_{r}$ is the temperature of the reservoir Seifert (2005); Qian (2001); Ge
and Qian (2010); Tomé and de Oliveira (2010). Together with the expression of
probability current (Eq.7), Eq.8 yields Seifert (2005)

$\dot{S}_{\text{\text{tot}}}(t)=\int\frac{\vec{J}^{\intercal}({\bf
x},t)\cdot{\bf D}^{-1}({\bf x})\cdot\vec{J}({\bf x},t)}{P({\bf x},t)}d{\bf
x}.$ (9)

The entropy production rate at steady state, obtained from $J^{\text{ss}}({\bf
x})$ and $P^{\text{ss}}({\bf x})$, allows us to evaluate the uncertainty
product $\mathcal{Q}$ at steady state (Eq.5).

Cautionary remarks are in place regarding the use of Eq.9 to evaluate the
entropy production, which is derived from the Fokker-Planck equation
describing the time evolution of the probability density for the dynamic
variables ${\bf x}$. In this study, we employed a widely adopted strategy of
approximating the chemical master equations (CME) for the Brusselator and the
PFK model for glycolytic oscillation to the corresponding Fokker-Planck
equations via the van Kampen’s linear noise approximation (LNA or
$\Omega$-expansion) (see SI) van Kampen (2007); Schuster (2016); Qian  _et
al._ (2002); Wang  _et al._ (2008); Xiao  _et al._ (2008); Cao  _et al._
(2015). However, consistency between CME and LNA approach has recently been
questioned in the context of stochastic thermodynamics Grima (2010); Grima
_et al._ (2011); Horowitz (2015). If the system size parameter $\Omega$ is too
small, not only the accuracy of the approximation becomes questionable Grima
(2010); Grima  _et al._ (2011), but the entropy production rate calculated by
employing Eq.9 is also bound to _underestimate_ the true value Horowitz
(2015). Since the system size we have chosen for the simulation
($\Omega=1600$) is large enough that the approximation is essentially taken in
the regime where the discrepancy between the entropy productions calculated
from CME and from Fokker-Planck approach should not be significant. Another
possible cause of underestimation of entropy production arises when one adopts
the coarse-graining Yu  _et al._ (2021) or the projection of dynamics to slow
degrees of freedom Zwanzig (2001); Van den Broeck and Esposito (2010). For the
cases of the Brusselator and the glycolytic oscillation studied here, their
reaction dynamics are defined with two stochastic variables ($x$ and $y$ for
the Brusselator; $[S]$ and $[P]$ for the glycolytic oscillation), we probe
both stochastic variables that are slowly varying with time and faithfully
represent the oscillatory dynamics.

Figure 1: A. Dynamic phase diagram of the Brusselator. The steady state
oscillation of limit cycle occurs in the pink region, left side of the phase
boundary. The blue region marked with “stable fixed points” displays no
oscillations at steady states. Two exemplary trajectories, one (red star) from
the region of limit cycles and the other (cyan star) from the region of stable
fixed points are shown. The corresponding phase planes with the vector fields
are depicted in Fig. S2. B. The 2D diagram of $\mathcal{Q}(a,b)$ corresponding
to the phase region of limit cycles. C. The 2D diagrams for other quantities,
$\dot{S}_{\text{tot}}$, $\langle T\rangle$, $\langle(\delta T)^{2}\rangle$,
$J_{\text{cycle}}$, and $A$, as a function of $a$ and $b$.

###  I.2 Brusselator

The Brusselator, a model for autocatalytic reactions that display sustained
oscillatory dynamics in certain range of parameters, offers all the
ingredients of oscillatory dynamics required to learn the dissipation and
precision of glycolytic oscillations to be studied.

The model forms an open system, consisting of two chemical compounds $X$ and
$Y$ reacting each other and being depleted while other source compounds $A$
and $B$ are constantly supplied to the system at fixed concentrations. The
reaction scheme for Brusselator is

$\displaystyle A$ $\displaystyle\xrightarrow{k_{1}}X,$ $\displaystyle B$
$\displaystyle\xrightarrow{k_{2}}Y,$ $\displaystyle 2X+Y$
$\displaystyle\xrightarrow{k_{3}}3X,$ $\displaystyle X$
$\displaystyle\xrightarrow{k_{4}}\phi.$ (10)

Properly non-dimensionalized (see SI), the rate equations at the limit of an
infinite volume $(\Omega\rightarrow\infty)$ is given as

$\displaystyle\frac{dx}{dt}$ $\displaystyle=f(x,y)=a-x+x^{2}y,$
$\displaystyle\frac{dy}{dt}$ $\displaystyle=g(x,y)=b-x^{2}y,$ (11)

where $x$, $y$, $a$, and $b$ are the non-dimensionalized concentrations of
chemical species $X$, $Y$, $A$ and $B$, respectively. The equations for the
concentrations of $x$ and $y$ are nonlinear; thus when $x$ and $y$ are
expanded around a fixed point $(x^{\ast},y^{\ast})$ that satisfies
$f(x^{\ast},y^{\ast})=0$ and $g(x^{\ast},y^{\ast})=0$ such that
$x=x^{\ast}+\delta x$ and $y=y^{\ast}+\delta y$, it yields a set of linearized
equations

$\displaystyle\delta\dot{{\bf
x}}=\mathcal{J}(x^{\ast},y^{\ast})\cdot\delta{\bf x}.$ (12)

where $\delta{\bf x}\equiv(\delta x,\delta y)^{\intercal}$ and

$\displaystyle\mathcal{J}(x^{\ast},y^{\ast})$
$\displaystyle=\begin{bmatrix}f_{x}&f_{y}\\\ g_{x}&g_{y}\end{bmatrix}$

is the Jacobian matrix evaluated at $(x^{\ast},y^{\ast})=(a+b,b/(a+b)^{2})$,
where $f_{x}\equiv\partial_{x}f(x,y)|_{(x,y)=(x^{\ast},y^{\ast})}$. Since the
fluctuation $\delta{\bf x}$ is expected to change with time as $\delta{\bf
x}\sim e^{\lambda t}$, the stability of the fixed point is determined by the
real parts of eigenvalues ($\lambda$) obtained from the characteristic
equation
$\det{\left[\lambda\mathcal{I}-\mathcal{J}(x^{\ast},y^{\ast})\right]}=0$ where
$\mathcal{I}$ is the identity matrix. The eigenvalues are:

$\displaystyle\lambda=\frac{1}{2}(\tau\pm\sqrt{\tau^{2}-4\Delta})$ (13)

with $\tau=-(a+b)^{2}-(a-b)/(a+b)$ and $\Delta=(a+b)^{2}$. Since $\Delta>0$
for any $a$ and $b$, the sign of ${\bf Re}(\lambda)$ is determined entirely by
the sign of $\tau$. The condition of $\tau=0$, leading to
$b^{3}+3ab^{2}+(3a^{2}-1)b+a^{3}+a=0$, determines the phase boundary (Fig.1A).
The set of parameters $(a,b)$ belonging to the blue region of phase diagram
(Fig.1A) leads to $\tau<0$, then the fixed points are stable, and the time
trajectory of $(x(t),y(t))$ converges to $(x^{\ast},y^{\ast})$ (see the lower
rightmost panel of Fig.1A). On the other hand, the region colored in red
($\tau>0$) yields unstable fixed points that produce limit cycles (see Fig.
S2A and B for the vector fields $(\dot{x}(t),\dot{y}(t))$ depicted for $(a,b)$
leading to unstable and stable fixed points). For the case of 2D phase plane,
existence of limit cycles is always guaranteed by the Poincaré-Bendixson
theorem Strogatz (2014).

Figure 2: Allosteric model for PFK in glycolysis. The PFK1 is an oligomeric
enzyme, which takes different oligomeric state depending on the organism.
Depicted is the octameric structure (dimer of tetramer) of PFK1 in $R$ state
(PDB code 4U1R Kloos  _et al._ (2015)). Each subunit has substrate ($S=$ F6P)
and product ($P=$ FBP) binding site, the latter of which allosterically
regulate the enzymatic activity of the enzyme. The PFK model of glycolytic
oscillation assumes an open thermodynamic system where glucose (or substrate)
molecules are injected at the rate $\nu$ and the product ($P$) is drained at
$k_{s}$. Depicted is the octameric form of enzyme with $n=8$ binding sites for
substrate and product. At low substrate concentration, the enzyme is in the
_tense_ ($T$) state. But, with increasing substrate concentration, transition
to the $R$ state occurs, increasing the catalytic activity of enzyme.

For a system with a finite volume, $\Omega$, the FP equation derived from CME
Qian  _et al._ (2002) for the Brusselator is obtained with the driving force
vector $\vec{F}(x,y)$ and diffusion tensor ${\bf D}(x,y)$.

$\displaystyle\vec{F}(x,y)$ $\displaystyle=\begin{bmatrix}a-x+x^{2}y\\\
b-x^{2}y\end{bmatrix}+\frac{1}{2\Omega}\begin{bmatrix}-1/2-2xy+x^{2}/2&\\\
2xy-x^{2}/2\end{bmatrix},$ $\displaystyle{\bf D}(x,y)$
$\displaystyle=\frac{1}{2\Omega}\begin{bmatrix}a+x+x^{2}y&-x^{2}y\\\
-x^{2}y&b+x^{2}y\end{bmatrix}.$ (14)

The $(a,b)$-dependent steady state probability distribution
$P_{\text{ss}}(x,y)$, obtained from an ensemble of trajectories generated
using the simulations based on Gillespie’s algorithm, allows us to calculate
the $\vec{J}_{\text{ss}}$ (Eq.6), and hence $\dot{S}_{\text{\text{tot}}}$ at
the steady state using Eq.9.

In the parameter range of $(a,b)$ yielding unstable fixed points ($\tau>0$),
we have computed $(a,b)$-dependent 2D diagrams of various quantities at
$\Omega=1600$: $\dot{S}_{\text{tot}}$, $\langle T\rangle$, $\langle\delta
T^{2}\rangle$, amplitude of oscillations ($A$), and the integral current
$J_{\text{cycle}}=\oint\vec{J}_{\text{ss}}\cdot d\vec{l}\Big{/}\oint dl$ (Fig
1). The entropy production rate $\dot{S}_{\text{tot}}$ show overall positive
correlation with $\langle T\rangle$, $\langle\delta T^{2}\rangle$, and $A$,
but not with $J_{\text{cycle}}$. The larger $J_{\text{cycle}}$ signifies
faster oscillations. It is noteworthy that the correlation or anti-correlation
between the quantities calculated here is relatively clear for the
Brusselator, but same is not necessarily true for glycolytic oscillator
(compare Fig. S3A and B). The $\langle\delta T^{2}\rangle$, displaying a non-
monotonic variation, is minimized at a basin of parameter space around
$(a,b)\approx(0.13,0.55)$. The product of $\dot{S}_{\text{tot}}$ and
$\langle\delta T^{2}\rangle/\langle T\rangle$ gives rise to the 2D diagram of
$\mathcal{Q}(a,b)$ (Fig.1), indicating that $\mathcal{Q}$ is minimized in the
vicinity of the phase boundary $(a,b)=(0.14,0.5)$ to $\mathcal{Q}\approx 17$.

Importantly, along a noisy limit cycle over the range of $\Omega$ being varied
(see Fig. S4), $\mathcal{Q}$ is independent of the system size $\Omega$. The
entropy production is an extensive quantity that linearly increase with
$\Omega$ for a given time interval $t$. Thus, the entropy production rate
scales with the volume as $\dot{S}_{\text{tot}}\sim\Omega$.Xiao  _et al._
(2008, 2009) Next, the fluctuation of the oscillatory period $\langle(\delta
T)^{2}\rangle$ is proportional to the magnitude of the $\Omega$-dependent
diffusion tensor ${\bf D}$, such that $\langle(\delta T)^{2}\rangle\sim{\bf
D}\sim\Omega^{-1}$ as defined in Eq.14, whereas $\langle T\rangle$ is decided
independently from $\Omega$. Taken together, the uncertainty product
$\mathcal{Q}$ is a quantity independent of $\Omega$, which can also be
confirmed using a host of simulations carried out at fixed parameter values
with varying $\Omega$ (see Fig. S4).

##  II Results

The allosteric regulation of PFK1 and its substrate and product concentration-
dependent enzymatic activity can be formulated using the strategy of Mono-
Wyman-Changeaux model Mono  _et al._ (1965); Thirumalai  _et al._ (2019). The
enzyme PFK1, a oligomeric complex consisting of $n$ catalytic and $n$
regulatory sites to which the substrate and product, respectively, can bind,
is equilibrated between tense (inactive, $T$) and relaxed (active, $R$)
states, which differ in terms of their conformations and binding affinities to
substrate and product (see Fig.2). $T$ state can only accommodate substrate
molecules, and the subscript $i(=0,1,\ldots n)$ in $T_{i}$ denotes the number
of substrates bound to the binding sites. On the other hand, $R$ state can
accommodate both substrate and product molecules; the two subscripts
$i(=0,1,\ldots,8)$ and $j(=0,1,\ldots,n)$ of $R_{ij}$ denote the number of
substrates and products bound to the catalytic and regulatory sites,
respectively. In the absence of substrate, two apo states of PFK1,
$R_{0}(\equiv R_{00})$ and $T_{0}$ states, are chemically equilibrated with
the ratio, $L=[T_{0}]/[R_{0}]$ called allosteric constant, which determines
the degree of cooperativity of the enzyme. The glucose converting into
substate $S$ (fructose-6-phosphate, F6P) via multiple steps is injected at a
constant rate $\nu$ while the product $P$ (fructose-1,6-biopohosphate, FBP)
either binds exclusively to the allosteric sites of the $R$ state acting as a
positive activator, or degrades at a rate of $k_{s}$. The increase of F6P as a
result of the catalytic processes of glycolysis in turn increases the amount
of FBP, which positively regulates the catalytic activity of PFK1 by promoting
the $T$-to-$R$ transitions. The nonlinear response of FBP-binding induced
autocatalytic activation of PFK1 generates the glycolytic oscillation Tornheim
(1988); Yaney  _et al._ (1995).

We assume that the binding and unbinding of substrate ($S$) and product ($P$)
to and from each binding site of $R$ state occur with the rates
$k_{R,\text{on}}^{S}$, $k_{R,\text{off}}^{S}$, and $k_{R,\text{on}}^{P}$,
$k_{R,\text{off}}^{P}$, and that only the $S$ can bind/unbind to $T$ state
with $k_{T,\text{on}}^{S}$, $k_{T,\text{off}}^{S}$. Then, the concentration of
each state of PFK1 is obtained as follows by assuming a quasi-steady-state
approximation Goldbeter and Lefever (1972).

$\displaystyle[T_{i}]$
$\displaystyle=\binom{n}{i}\left(\frac{[S]}{K_{T}^{S}}\right)^{i}[T_{0}]=L\binom{n}{i}\left(\frac{c[S]}{K_{R}^{S}}\right)^{i}[R_{0}],$
$\displaystyle[R_{ij}]$
$\displaystyle=\binom{n}{i}\left[\alpha([S])\right]^{i}\binom{n}{j}\left[\gamma([P])\right]^{j}[R_{0}]$
(15)

where

$\displaystyle\alpha([S])$
$\displaystyle=\frac{[S]}{(K_{R}^{S}+k/k_{R,\text{on}}^{S})},$
$\displaystyle\gamma([P])$ $\displaystyle=\frac{[P]}{K_{R}^{P}},$

$K_{R}^{S}(=k_{R,\text{off}}^{S}/k_{R,\text{on}}^{S})$ and
$K_{T}^{S}(=k_{T,\text{off}}^{S}/k_{T,\text{on}}^{S})$ are the binding
affinities (dissociation constants) of the substrate to the catalytic site in
the $R$ and $T$ states, respectively, whereas $K_{R}^{P}$ is the binding
affinity of the product to the regulatory site in the $R$ state. The parameter
$c=K_{R}^{S}/K_{T}^{S}$ denotes the ratio of the dissociation constants of the
substrate from the catalytic sites in $R$ and $T$ states. Then, from the
concentration of each enzyme state $[T_{i}]$ and $[R_{ij}]$ (Eq.15), it is
straightforward to calculate a binding polynomial $\bar{Y}$ ($0\leq\bar{Y}\leq
1$), namely the fraction of catalytic sites in the $R$ state bound by the
substrate,

$\displaystyle\bar{Y}$ $\displaystyle=\frac{\text{total substrates bound to
$R$ state}}{n\times\text{total enzymes}}$
$\displaystyle=\frac{\sum_{i=0}^{n}\sum_{j=0}^{n}i[R_{ij}]}{nZ([S],[P])}$
$\displaystyle=\frac{\alpha([S])q([S],[P])}{nZ([S],[P])}$ (16)

where

$\displaystyle q([S],[P])$ $\displaystyle\equiv
n(1+\alpha([S]))^{(n-1)}(1+\gamma([P]))^{n}[R_{0}],$

and the total concentration of enzyme $Z([S],[P])$

$\displaystyle Z$
$\displaystyle([S],[P])=\sum_{i=0}^{n}[T_{i}]+\sum_{i=0}^{n}\sum_{j=0}^{n}[R_{ij}]$
$\displaystyle=\left[L\left(1+\frac{c[S]}{K_{R}^{S}}\right)^{n}+(1+\alpha([S]))^{n}(1+\gamma([P]))^{n}\right][R_{0}]$
(17)

Figure 3:  A. Phase diagram of glycolytic oscillations as a function of $\nu$
and $k_{s}$ with the allosteric constant ($L=4\times 10^{9}$), catalytic rate
($k=500$ $s^{-1}$) and $c=0.01$. No catalysis occurs in the hashed region. The
simulations were performed with $\nu=0.005-0.03$ mM/s, demarcated with the
dark pink region, corresponding to the glucose uptake rate of yeast. B. Three
representative trajectories showing oscillatory dynamics of substrate and
product concentrations. The trajectory marked with yellow star was generated
at the parameter values relevant for octameric PFK of yeast giving rise to the
oscillatory period of $\sim$ 400 s Blangy  _et al._ (1968); Hess  _et al._
(1969); Goldbeter and Lefever (1972); McCoy  _et al._ (2005); Cardon and Boyer
(1978); Moyer  _et al._ (1998) ($L=4\times 10^{9}$, $k=500$ s-1, $\nu=0.005$
mM/s, $k_{s}=0.05$ s-1, $c=0.01$). C. 2D diagram of uncertainty product,
$\mathcal{Q}$ as a function of $\nu$ and $k_{s}$. D. $\dot{S}_{\text{tot}}$,
$\langle T\rangle$, $\langle(\delta T)^{2}\rangle$, $J_{\text{cycle}}$, and
$A$ calculated as a function of $\nu$ and $k_{s}$.

Substrates supplied with constant rate $\nu$ bind to $T$ or $R$ state, and is
catalyzed by the $R$ state of PFK1 ($R_{ij}$) with a rate $k$ to generate
product molecules, whereas the products either bind to the regulatory site of
$R$ state or depleted from the system at a rate $k_{s}$. The mass action laws
of the substrate and product yield a set of coupled nonlinear equations
Goldbeter and Lefever (1972):

$\displaystyle\frac{d[S]}{dt}$
$\displaystyle=\nu-k\sum_{i=0}^{n}\sum_{j=0}^{n}i[R_{ij}]=\nu-k\alpha([S])q([S],[P]),$
$\displaystyle\frac{d[P]}{dt}$
$\displaystyle=k\sum_{i=0}^{n}\sum_{j=0}^{n}i[R_{ij}]-k_{s}[P]=k\alpha([S])q([S],[P])-k_{s}[P].$
(18)

Here, by specifically considering the yeast PFK which adopts an octameric form
($n=8$), we use the allosteric constant $L=4\times 10^{9}$ which is 3 orders
of magnitude greater than the value known for the tetramer ($n=4$) Blangy  _et
al._ (1968). Other parameters related with substrate/product binding to each
protomer are expected to be identical between tetramer and octamer, thus we
take $k_{R,\text{on}}^{S}=2000$/mM$\cdot$s, $K^{S}_{R}=0.05$ mM and
$K^{P}_{R}=0.025$ mM McCoy  _et al._ (2005); Moyer  _et al._ (1998); Blangy
_et al._ (1968); Cardon and Boyer (1978). The ratio of the substrate binding
affinities to the $R$ and $T$ states is $c=K_{R}/K_{T}=0.01$, so that the
substrate binds preferentially to the $R$ state. The ATP hydrolysis time due
to ATPase activity is typically $\gtrsim\mathcal{O}(1)$ msec Gilbert and
Johnson (1994), and hence we set $k=500$ s-1.

Following the same procedure used in the analysis of Brusselator, we analyze
Eq.18 to calculate phase diagrams as a function of $\nu$ and $k_{s}$ (Fig.3A)
and generate dynamical trajectories of $[S](t)$ and $[P](t)$ (Fig.3B). For the
values of $(\nu,k_{s})$ pertaining to the limit cycles, the substrate
concentration display saw-tooth like oscillatory pattern in time, and and the
product concentration spikes when the substrate concentration falls, which
generates a loop in the phase plane of $([S],[P])$ (Fig.3B). A trajectory with
the mean oscillatory period of $\sim$ 400 s ($\approx 6-7$ min), fluctuating
between $\sim$ 0.1 mM and $\sim$ 1.5 mM, emerges at the condition ($L=4\times
10^{9}$, $\nu=0.005$ mM/s, $k_{s}=0.05$ s-1, $c=0.01$) (yellow stars in Figs.
3A, B, C). The period and amplitude of the oscillation comport well with those
observed in yeast and yeast extract Hess  _et al._ (1969); Hess and Szabo
(1979), in which PFK1 enzymes are oligomerized to an octameric form.

Next, the FP equation for the glycolytic oscillations is obtained from the CME
corresponding to Eq.18 with the following force vector and diffusion tensor:

$\displaystyle\vec{F}$ $\displaystyle=\begin{bmatrix}\nu-k\alpha q\\\ k\alpha
q-k_{s}[P]\end{bmatrix}+\frac{1}{2\Omega}\begin{bmatrix}-k\alpha
q-k\alpha\partial_{[S]}q+k\alpha\partial_{[P]}q&\\\ -k_{s}+k\alpha
q+k\alpha\partial_{[S]}q-k\alpha\partial_{[P]}q\end{bmatrix},$
$\displaystyle{\bf D}$
$\displaystyle=\frac{1}{2\Omega}\begin{bmatrix}\nu+k\alpha q&-k\alpha q\\\
-k\alpha q&k_{s}[P]+k\alpha q\end{bmatrix},$ (19)

where the concentration dependences of $\alpha=\alpha([S])$ (Eq.II) and
$q=q([S],[P])$ (Eq.II) are omitted for the simplicity of the expression. The
FP equation with these $\vec{F}$ and ${\bf D}$ is used to calculate
$J^{\text{ss}}([S],[P])$, $P^{\text{ss}}([S],[P])$, and $\dot{S}_{\text{tot}}$
based on Eq.9.

Shown in Fig.3C, D are the 2D diagrams of $\dot{S}_{\text{tot}}$, $\langle
T\rangle$, $\langle\delta T^{2}\rangle$, $J_{\text{cycle}}$, $A$, and finally
$\mathcal{Q}$ as a function of $\nu$ and $k_{s}$, which are the two
experimentally controllable parameters. Overall, the correlations between
these quantities are not so strong in comparison with those calculated for the
Brusselator (see Fig. S3). It is fair to say that the correlation or a trend
seen in Brusselator cannot be generalized to other biochemical oscillators. At
the parameter values, $\nu=0.005$ mM/s-1 and $k_{s}=0.05$ s-1, yielding
$T=\langle T\rangle\pm\langle\delta T^{2}\rangle^{1/2}\approx 400\pm 20$ s,
the uncertainty product is $\mathcal{Q}\simeq 31$. Remarkably, while
$\mathcal{Q}\simeq 31$ is observed in the vicinity of the lower phase boundary
where the entropy production rate is minimal over the relevant phase space
(see Fig.3C and the first panel of Fig.3D), $\langle T\rangle\approx 400$ sec
(the second panel of Fig.3D) is obtained only at
$(\nu,k_{s})\simeq(0.005\text{ mM/s},0.01\text{ s}^{-1})$.

##  III Discussions

Glycolytic oscillations are the temporal order that emerges under certain
special conditions in which parameters defining a set of coupled nonlinear
equations yield unstable fixed points. Notably, for octameric form of PFK
oligomers to demonstrate oscillatory dynamics, the condition of
$c=K_{R}/K_{T}=0.01$, which renders the substrate binding to the protomer in
$R$ state more preferable than to $T$ state by a hundred fold, is essential.
If $c$ is increased to 0.1, the phase space corresponding to limit cycles
significantly narrows down (Fig.4). For tetrameric form of PFK ($n=4$), which
pertains to bacteria, the phase space region for limit cycles is much narrower
even when $c=0.01$ and $L$ is set to the value of octamer ($L=4\times 10^{9}$)
(Fig.4). Unless $k_{s}$ is tuned to a narrow interval of $k_{s}\simeq
0.01-0.1$ s-1, no oscillation is expected. The $(\nu,k_{s})$ phase diagrams of
glycolysis with varying $c$ and $n$ (Fig.4) rationalize why glycolytic
oscillations were only reported in eukaryotic PFK, where PFK exists in the
octameric form.

The TUR, which specifies the physical lower bound to the uncertainty product,
is used to assess how the period of temporal order emerging from the
underlying dynamical process is balanced with the dissipation under the
constraint of cost-precision trade-off. In the Brusselator, whose dynamical
behavior is defined only with two parameters ($a$ and $b$), there is a
specific case that both precision of oscillatory period and dissipation are
simultaneously minimized to yield a reasonably small uncertainty product
$\mathcal{Q}\simeq 17$ over the phase space. In comparison, glycolytic
oscillations are more complicated with many more parameters ($L$, $c$, $n$,
$K_{R}^{S}$, $K_{R}^{P}$, $k$, $\nu$, $k_{s}$). To simplify the problem, we
have reduced the unknowns by assuming that some of the parameters have
identical values with those pertaining to the protomer. The values of
uncertainty product $\mathcal{Q}$ for the glycolytic oscillations at their
working condition producing the period of $\sim$ (5 – 10) min is
$\mathcal{Q}\simeq 31$. Remarkably, given the substrate injection rate
$\nu=0.005-0.05$ s-1, $\mathcal{Q}\simeq 31$ is effectively the minimal value
over the phase space involving the limit cycle (Fig.3C). $\mathcal{Q}\simeq
31$ is greater than those determined for the molecular motors
$\mathcal{Q}\approx 7-15$ Hwang and Hyeon (2018); Mugnai  _et al._ (2020), and
biological copy process by exonuclease-deficient T7 DNA polymerase
$\mathcal{Q}\approx 10$ Song and Hyeon (2020), but smaller than
$\mathcal{Q}\approx 45-50$ for the translation process by E. coli ribosome
Piñeros and Tlusty (2020); Song and Hyeon (2020). In comparison with the
uncertainty product determined for other biochemical cycles
($\mathcal{Q}\approx 10^{3}$) Marsland III  _et al._ (2019), which severely
underperform the TUR’s lower bound of $2$, the value of the uncertainty
product $\mathcal{Q}\simeq 31$ for the glycolytic oscillation arising from
octameric PFK is significantly smaller, minimizing the entropy production rate
over the relevant phase, which indicates the cost-effectiveness of the
molecular mechanism generating the oscillatory dynamics.

Figure 4: The 2D phase diagram of glycolytic oscillations as a function of the
injection ($\nu$) and degradation rates ($k_{s}$) for varying $c$ at
$n(=2,4,8)$ with $L=4\times 10^{9}$ and $k=500$ s-1. No catalysis occurs in
the hashed region. The working condition of glycolytic oscillations
($\nu=0.005$ mM/s, $k_{s}=0.05$ s-1) that produces the oscillations with the
mean period $\langle T\rangle=400$ s at $c=0.01$ and $n=8$ (Fig.3) is marked
with the yellow star in each panel.

Lastly, it is of particular note that the $\sim$ (5 – 10) min oscillation
period is observed in cellular or physiological scales as well, such as the
blood glucose level, intracellular Ca2+ concentrations, and membrane action
potentials, which are controlled by the pulsatile secretion of insulin with a
period of $\sim$ (5 – 10) min McKenna  _et al._ (2016); Lang  _et al._ (1979);
Westermark and Lansner (2003), suggestive of a connection between the dynamics
at the molecular and macroscopic scales Bertram  _et al._ (2004), and their
synchronization Lee  _et al._ (2017). Remarkably, the loss of pulsatile
insulin release resulting from elevated glucose level McKenna  _et al._
(2016); Lee  _et al._ (2017) is also consistent with our study that the
working condition of glycolytic oscillation is situated near the borderline of
the phase boundary. The phase diagram depicted in Fig.3A predicts that a
moderate elevation of glucose injection rate beyond $\nu\approx 0.01$ mM/s for
a fixed $k_{s}=0.05$ s-1 would abolish the oscillations.

###### Acknowledgements.

We thank Prof. Junghyo Jo for illuminating discussions on glucose level
oscillations. This study is supported by KIAS Individual Grants CG076501
(P.K.) and CG035003 (C.H.). We thank the Center for Advanced Computation in
KIAS for providing computing resources.

## References

  * Bustamante  _et al._ (2005) C. Bustamante, J. Liphardt, and F. Ritort, Physics Today 58, 43 (2005).
  * Rust  _et al._ (2011) M. J. Rust, S. S. Golden, and E. K. O’Shea, Science 331, 220 (2011).
  * Rust  _et al._ (2007) M. J. Rust, J. S. Markson, W. S. Lane, D. S. Fisher, and E. K. O’Shea, Science 318, 809 (2007).
  * Novak and Tyson (1993) B. Novak and J. Tyson, J. Cell Sci. 106, 1153 (1993).
  * Tyson and Novak (2001) J. Tyson and B. Novak, J. Theor. Biol. 210, 249 (2001).
  * Ferrell  _et al._ (2011) J. E. Ferrell, T. Y.-C. Tsai, and Q. Yang, Cell 144, 874 (2011).
  * Prigogine (1978) I. Prigogine, Science 201, 777 (1978).
  * Goldbeter (2018) A. Goldbeter, Phil. Trans. Roy. Soc. A: Math. Phys. Eng. Sci. 376, 20170376 (2018).
  * Duysens and Amesz (1957) L. Duysens and J. Amesz, Biochimi. Biophys. Acta 24, 19 (1957).
  * Boiteux  _et al._ (1975) A. Boiteux, A. Goldbeter, and B. Hess, Proc. Natl. Acad. Sci. U. S. A. 72, 3829 (1975).
  * Goldbeter and Caplan (1976) A. Goldbeter and S. R. Caplan, Annu. Rev. Biophys. Bioeng. 5, 449 (1976).
  * Tornheim (1997) K. Tornheim, Diabetes 46, 1375 (1997).
  * Bertram  _et al._ (2004) R. Bertram, L. Satin, M. Zhang, P. Smolen, and A. Sherman, Biophys. J. 87, 3074 (2004).
  * Hess and Boiteux (1971) B. Hess and A. Boiteux, Annu. Rev. Biochem. 40, 237 (1971).
  * Barato and Seifert (2015) A. C. Barato and U. Seifert, Phys. Rev. Lett. 114, 158101 (2015).
  * Gingrich  _et al._ (2016) T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Phys. Rev. Lett. 116, 120601 (2016).
  * Horowitz and Gingrich (2019) J. M. Horowitz and T. R. Gingrich, Nat. Phys. 16, 15 (2019).
  * Hasegawa and Van Vu (2019) Y. Hasegawa and T. Van Vu, Phys. Rev. Lett. 123, 110602 (2019).
  * Proesmans and Van den Broeck (2017) K. Proesmans and C. Van den Broeck, E 119, 20001 (2017).
  * Lee  _et al._ (2019) J. S. Lee, J.-M. Park, and H. Park, Phys. Rev. E 100, 062132 (2019).
  * Agarwalla and Segal (2018) B. K. Agarwalla and D. Segal, Phys. Rev. B 98, 155438 (2018).
  * Lee  _et al._ (2018) S. Lee, C. Hyeon, and J. Jo, Phys. Rev. E 98, 032119 (2018).
  * Potts and Samuelsson (2019) P. P. Potts and P. Samuelsson, Phys. Rev. E 100, 052137 (2019).
  * Koyuk  _et al._ (2018) T. Koyuk, U. Seifert, and P. Pietzonka, J. Phys. A: Math. Theor. 52, 02LT02 (2018).
  * Koyuk and Seifert (2019) T. Koyuk and U. Seifert, Phys. Rev. Lett. 122, 230601 (2019).
  * Koyuk and Seifert (2020) T. Koyuk and U. Seifert, Phys. Rev. Lett. 125, 260604 (2020).
  * Pietzonka  _et al._ (2016) P. Pietzonka, A. C. Barato, and U. Seifert, J. Stat. Mech. Theory Exp. , 124004 (2016).
  * Hwang and Hyeon (2018) W. Hwang and C. Hyeon, J. Phys. Chem. Lett. 9, 513 (2018).
  * Uhl and Seifert (2018) M. Uhl and U. Seifert, Phys. Rev. E 98, 022402 (2018).
  * Song and Hyeon (2020) Y. Song and C. Hyeon, J. Phys. Chem. Lett. 11, 3136 (2020).
  * Piñeros and Tlusty (2020) W. D. Piñeros and T. Tlusty, Phys. Rev. E 101, 022415 (2020).
  * Marsland III  _et al._ (2019) R. Marsland III, W. Cui, and J. M. Horowitz, J. Roy. Soc. Interface 16, 20190098 (2019).
  * Song and Hyeon (2021) Y. Song and C. Hyeon, J. Chem. Phys. 154, 130901 (2021).
  * Gingrich and Horowitz (2017) T. R. Gingrich and J. M. Horowitz, Phys. Rev. Lett. 119, 170601 (2017).
  * Schnitzer and Block (1995) M. J. Schnitzer and S. Block, Cold spring harbor symposia on quantitative biology 60, 793 (1995).
  * Cao  _et al._ (2015) Y. Cao, H. Wang, Q. Ouyang, and Y. Tu, Nat. Phys. 11, 772 (2015).
  * Morelli and Jülicher (2007) L. G. Morelli and F. Jülicher, Phys. Rev. Lett. 98, 228101 (2007).
  * Marsland III and England (2017) R. Marsland III and J. England, Rep. Prog. Phys. 81, 016601 (2017).
  * Seifert (2005) U. Seifert, Phys. Rev. Lett. 95, 040602 (2005).
  * Qian (2001) H. Qian, Phys. Rev. E 64, 022101 (2001).
  * Ge and Qian (2010) H. Ge and H. Qian, Phys. Rev. E 81, 051133 (2010).
  * Tomé and de Oliveira (2010) T. Tomé and M. J. de Oliveira, Phys. Rev. E 82, 021120 (2010).
  * van Kampen (2007) N. G. van Kampen,  _Stochastic Processes in Chemistry and Physics_ (Elsevier, North Holland, Amsterdam, 2007).
  * Schuster (2016) P. Schuster,  _Stochasticity in processes_ (Springer, Berlin, 2016).
  * Qian  _et al._ (2002) H. Qian, S. Saffarian, and E. L. Elson, Proc. Natl. Acad. Sci. U. S. A. 99, 10376 (2002).
  * Wang  _et al._ (2008) J. Wang, L. Xu, and E. Wang, Proc. Natl. Acad. Sci. U. S. A. 105, 12271 (2008).
  * Xiao  _et al._ (2008) T. J. Xiao, Z. Hou, and H. Xin, J. Chem. Phys. 129, 114506 (2008).
  * Grima (2010) R. Grima, J. Chem. Phys. 133, 07B604 (2010).
  * Grima  _et al._ (2011) R. Grima, P. Thomas, and A. V. Straube, J. Chem. Phys. 135, 084103 (2011).
  * Horowitz (2015) J. M. Horowitz, J. Chem. Phys. 143, 044111 (2015).
  * Yu  _et al._ (2021) Q. Yu, D. Zhang, and Y. Tu, Phys. Rev. Lett. 126, 080601 (2021).
  * Zwanzig (2001) R. Zwanzig,  _Nonequilibrium Statistical Mechanics_ (Oxford University press, New York, 2001).
  * Van den Broeck and Esposito (2010) C. Van den Broeck and M. Esposito, Phys. Rev. E 82, 011144 (2010).
  * Strogatz (2014) S. H. Strogatz,  _Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering_ (Westview Press, Boulder, 2014).
  * Xiao  _et al._ (2009) T. Xiao, Z. Hou, and H. Xin, J. Phys. Chem. B 113, 9316 (2009).
  * Mono  _et al._ (1965) J. Mono, J. Wyman, and J. P. Changeux, J. Mol. Biol. 12, 88 (1965).
  * Thirumalai  _et al._ (2019) D. Thirumalai, C. Hyeon, P. I. Zhuravlev, and G. H. Lorimer, Chem. Rev. 119, 6788 (2019).
  * Tornheim (1988) K. Tornheim, J. Biol. Chem. 263, 2619 (1988).
  * Yaney  _et al._ (1995) G. C. Yaney, V. Schultz, B. A. Cunningham, G. A. Dunaway, B. E. Corkey, and K. Tornheim, Diabetes 44, 1285 (1995).
  * Goldbeter and Lefever (1972) A. Goldbeter and R. Lefever, Biophys. J. 12, 1302 (1972).
  * Blangy  _et al._ (1968) D. Blangy, H. Buc, and J. Monod, J. Mol. Biol. 31, 13 (1968).
  * McCoy  _et al._ (2005) M. A. McCoy, M. M. Senior, and D. F. Wyss, J. Am. Chem. Soc. 127, 149 (2005).
  * Moyer  _et al._ (1998) M. Moyer, S. Gilbert, and K. Johnson, Biochemistry 37, 800—813 (1998).
  * Cardon and Boyer (1978) J. W. Cardon and P. D. Boyer, Eur. J. Biochem. 92, 443 (1978).
  * Gilbert and Johnson (1994) S. P. Gilbert and K. A. Johnson, Biochemistry 33, 1951 (1994).
  * Hess  _et al._ (1969) B. Hess, A. Boiteux, and J. Krüger, Adv. Enzyme Reg. 7, 149 (1969).
  * Hess and Szabo (1979) V. Hess and A. Szabo, J. Chem. Edu. 56, 289 (1979).
  * Mugnai  _et al._ (2020) M. L. Mugnai, C. Hyeon, M. Hinczewski, and D. Thirumalai, Rev. Mod. Phys. 92, 025001 (2020).
  * McKenna  _et al._ (2016) J. P. McKenna, R. Dhumpa, N. Mukhitov, M. G. Roper, and R. Bertram, PLoS Comp. Biol. 12, e1005143 (2016).
  * Lang  _et al._ (1979) D. A. Lang, D. R. Matthews, J. Peto, and R. C. Turner, New Eng. J. Med. 301, 1023 (1979).
  * Westermark and Lansner (2003) P. O. Westermark and A. Lansner, Biophys. J. 85, 126 (2003).
  * Lee  _et al._ (2017) B. Lee, T. Song, K. Lee, J. Kim, S. Han, P.-O. Berggren, S. H. Ryu, and J. Jo, PLoS one 12, e0172901 (2017).
  * Kloos  _et al._ (2015) M. Kloos, A. Bruser, J. Kirchberger, T. Schoeneberg, and N. Strater, Biochem. J. 469, 421 (2015).

##  IV Supporting Information

Non-dimensionalization of rate equations.  The stochastic version of the
Brusselator is written as

$\displaystyle\frac{dX}{dt}$
$\displaystyle=k_{1}A-k_{4}X+k_{3}\frac{X(X-1)Y}{\Omega^{2}},$
$\displaystyle\frac{dY}{dt}$
$\displaystyle=k_{2}B-k_{3}\frac{X(X-1)Y}{\Omega^{2}}$ (S1)

where $\Omega$ is the volume of the system. Using the following
transformations of the variables and parameters,

$\displaystyle\left(\frac{k_{3}}{k_{4}}\right)^{1/2}\frac{X}{\Omega}$
$\displaystyle\longrightarrow x,$
$\displaystyle\left(\frac{k_{3}}{k_{4}}\right)^{1/2}\frac{Y}{\Omega}$
$\displaystyle\longrightarrow y,$ $\displaystyle k_{4}t$
$\displaystyle\longrightarrow t,$
$\displaystyle\frac{k_{1}}{k_{4}}\left(\frac{k_{3}}{k_{4}}\right)^{1/2}\frac{A}{\Omega}$
$\displaystyle\longrightarrow a,$
$\displaystyle\frac{k_{2}}{k_{4}}\left(\frac{k_{3}}{k_{4}}\right)^{1/2}\frac{B}{\Omega}$
$\displaystyle\longrightarrow b,$ (S2)

one can write down a non-dimensionalized version of the rate equations at the
limit of $\Omega\rightarrow\infty$ as in the main text (Eq.11).

Diffusion approximation: Fokker Planck equation from Chemical Master
Equations. For a chemical species (X) involved in a reaction:

$sX+\cdots{\xrightarrow{k}}s^{\prime}X+\cdots$ (S3)

the time evolution equation of $X$ at the deterministic limit can be written
as

$\frac{d[X]}{dt}=Sk[X]^{s}$ (S4)

where $[X]$ is the concentration of $X$, $S=s^{\prime}-s$ is the
stoichiometric coefficient. The chemical state of the system at any time is
fully determined by the state vector ${\bf X}=(X_{1},X_{2},\dots,X_{N})$ where
$X_{i}$ is the number of chemical species $X_{i}$ in a compartment of a finite
volume $\Omega$. Then, probability distribution for the system to be in state
X at time $t$ is

$\displaystyle P(\textbf{X},t+dt)=P(\textbf{X},t)$
$\displaystyle+dt\sum_{r=1}^{R}\left[f_{r}(\textbf{X}-\textbf{S}_{r})P(\textbf{X}-\textbf{S}_{r},t)-f_{r}(\textbf{X})P(\textbf{X},t)\right].$
(S5)

where
$f_{r}(\textbf{X})=k_{r}\Omega\prod_{i=1}^{N}\frac{X_{i}!}{(X_{i}-s_{ir})!\Omega^{s_{ir}}}$
is the probability for the reaction $r$ to occur. At the limit of
$dt\rightarrow 0$, Eq.S5 is cast into the Chemical Master Equation (CME)

${\partial_{t}}P(\textbf{X},t)=\sum_{r=1}^{R}\left\\{f_{r}(\textbf{X}-\textbf{S}_{r})P(\textbf{X}-\textbf{S}_{r},t)-f_{r}(\textbf{X})P(\textbf{X},t)\right\\}.$
(S6)

Since analytic solutions of CME is known only for limited cases, CME is
typically approximated to Fokker Planck equation through a Taylor expansion of
the relevant terms to the second order,

$\displaystyle f_{r}(\textbf{X}-\textbf{S}_{r})$ $\displaystyle
P(\textbf{X}-\textbf{S}_{r},t)-f_{r}(\textbf{X})P(\textbf{X},t)$
$\displaystyle\approx-\sum_{i=1}^{N}S_{ir}\partial_{X_{i}}[f_{r}(\textbf{X})P(\textbf{X},t)]$
$\displaystyle+\sum_{i,j=1}^{N}S_{ir}S_{jr}\partial_{X_{i}}\partial_{X_{j}}[f_{r}(\textbf{X})P(\textbf{X},t)],$
(S7)

leading to

$\displaystyle{\partial_{t}}P(\textbf{X},t)$
$\displaystyle=-\sum_{i=1}^{N}\partial_{X_{i}}[A_{i}(\textbf{X})P(\textbf{X},t)]$
$\displaystyle+\frac{1}{2}\sum_{i,j=1}^{N}\partial_{X_{i}}\partial_{X_{j}}[B_{ij}(\textbf{X})P(\textbf{X},t)]$
(S8)

where the drift vector $A$ and diffusion matrix $B$ are given by

$\displaystyle A_{i}(\textbf{X})$
$\displaystyle=\sum_{r=1}^{R}S_{ir}f_{r}(\textbf{X}),$ $\displaystyle
B_{ij}(\textbf{X})$
$\displaystyle=\sum_{r=1}^{R}S_{ir}S_{jr}f_{r}(\textbf{X}).$ (S9)

In order to handle the time evolution of different chemical species inside the
compartment of volume $\Omega$ in terms of their concentrations, we define
$x_{i}\equiv X_{i}/\Omega\equiv[X_{i}]$. With this definition, the probability
density of a set of concentrations can be converted to that of molecular
counts as $P(\textbf{x},t)\equiv\Omega^{N}P(\Omega\textbf{x},t)$. Finally, the
Fokker-Planck equation for the stochastic time evolution of concentrations of
chemical species is obtained as

$\displaystyle{\partial_{t}}P(\textbf{x},t)$
$\displaystyle=-\sum_{i=1}^{N}\partial_{x_{i}}[A_{i}(\textbf{x})P(\textbf{x},t)]$
$\displaystyle+\frac{1}{2\Omega}\sum_{i,j=1}^{N}\partial_{x_{i}}\partial_{x_{j}}[B_{ij}(\textbf{x})P(\textbf{x},t)].$
(S10)

Figure S1: A. $\langle\delta T_{i}\delta T_{j}\rangle=\langle
T_{i}T_{j}\rangle-\langle T_{i}\rangle\langle T_{j}\rangle$ versus $|i-j|$.
When $|i-j|\gtrsim 3$, the correlation between two distinct periods is
effectively zero ($\langle\delta T_{i}\delta T_{j}\rangle\approx 0$). B.
$\langle\delta t^{2}\rangle$ versus $n$. The plot of $\langle\delta
t^{2}\rangle$ versus $n$ indicates that the variance of mean first passage
time $\langle\delta t^{2}\rangle$ increases linearly with $n$ (the number of
oscillations), satisfying $\langle\delta t^{2}\rangle=Dn$, where
$D=\langle\delta T^{2}\rangle$.  Figure S2: Vector field $(dx/dt,dy/dt)$ of
Brusselator around the fixed point depicted on the phase plane. (A) Unstable
fixed point with $a=0.15$ and $b=0.5$. (B) Stable fixed point with $a=0.28$
and $b=0.8$. $y=1/x-a/x^{2}$ and $y=b/x^{2}$ are drawn in pink and pale blue,
respectively.  Figure S3: Correlations between $\dot{S}_{\text{tot}}$ and
$\langle T\rangle$, $\langle(\delta T)^{2}\rangle$, $J_{\text{cycle}}$, $A$
for (A) brusselator and (B) glycolytic oscillator models. The data points were
generated by randomly selecting the parameter values pertaining to the phase
region of limit cycle.  Figure S4: Effect of system size on the entropy
production, fluctuations, oscillatory period, and $\mathcal{Q}$. The data are
fitted to $\dot{S}_{\text{tot}}\sim\Omega$ , $\langle\delta
T^{2}\rangle\sim\Omega^{-1}$. $\langle T\rangle$ and $\mathcal{Q}$ are
effectively constant, independent of $\Omega$.