Thermodynamic optimality of glycolytic oscillations
TThermodynamic optimality of glycolytic oscillations
Pureun Kim and Changbong Hyeon ∗ Korea Institute for Advanced Study, Seoul 02455, Korea
Temporal order in living matters reflects the self-organizing nature of dynamical processes driven out ofthermodynamic equilibrium. Because of functional reason, the period of a biochemical oscillation must betuned 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 reviewingthe basics of chemical oscillations using Brusselator as a model system, we study the glycolytic oscillationgenerated by octameric phosphofructokinase (PFK), which is known to display a period of several minutes. Byexploring the phase space of glycolytic oscillations, we find that the glycolytic oscillation under the cellularcondition is realized in a cost effective manner. Specifically, over the biologically relevant range of parametervalues of glycolysis and octameric PFK, the entropy production from the glycolytic oscillation is minimal whenthe oscillation period is (5 – 10) minutes. Further, the glycolytic oscillation is found at work near the phaseboundary of limit cycles, suggesting that a moderate increase of glucose injection rate leads to the loss ofoscillatory dynamics, which is reminiscent of the loss of pulsatile insulin release resulting from elevated bloodglucose level.
From the perspective of thermodynamics on themacroscale, it is clear that living organisms are sustained farfrom thermodynamic equilibrium in that energy and materialcurrents are supplied to the systems and dissipated to theenvironment [1]. For human consisting of ∼ trillions ofcells, the energy consumption is ∼ W. At the cellularscale, this amounts to as much as ∼ k B T of energybeing consumed per second, which drives a vast number ofcellular processes. Among them, biochemical oscillations,exemplified with circadian rhythm [2, 3] and the cell cycle[4–6], are the temporal dissipative structures that emerge frommolecular components and their dynamics in nonequilibrium[7, 8].The full catabolism of glucose under aerobic condition isthe primary source of ATP; one glucose molecule could yieldas many as 30-32 ATP molecules. Besides its utmost impor-tance along the metabolic pathway, the sustained oscillationof substrate and product concentrations with time period of ∼ (5 – 10) min observed in the glycolysis of yeast, muscle, andpancreatic β cell is a phenomenon of great interest on its own[9–13]. Studies carried out on yeast indicate that among a hostof glycolytic enzymes, the allostery of phosphofructokinase(PFK), which uses ATP hydrolysis to catalyze the phospho-rylation of fructose-6-phosphate (substrate, S ) into fructose-1,6-biopohosphate (product, P ), was identified to lie at thecore of the endogenous oscillatory mechanism [10, 14]. Thebinding of product to the allosteric site induces a cooperativetransition of the PFK enzyme from its inactive to active state,introducing a nonlinear positive response to the biochemicalcircuit. The product-activated allosteric regulation of PFK en-zyme becomes the source of the glycolytic oscillation.In this paper, we study the phosphofructokinase (PFK)model for glycolytic oscillations, and discuss how optimal thisenergy consuming process is in achieving its desired oscil-lation period at the expense of thermodynamic cost, whichwe quantify using the thermodynamic uncertainty relation(TUR), a recently discovered thermodynamic principle [15–18]. In the section Theoretical background , we briefly review the TUR for dynamical processes displaying temporal oscil-lations. Before studying the PFK model of glycolytic oscilla-tion, we review the Brusselator, a model chemical processesexhibiting oscillations, and analyze it in detail at varying pa-rameter values by means of the phase diagram to identify theoptimal condition of the process in light of TUR. The
Re-sults section covers the allosteric model of octameric PFK,dynamic phase diagram, and the mass balance equations ofPFK-catalyzed substrate and product giving rise to the gly-colytic oscillations. We assess the optimality of glycolytic os-cillations realized by octameric PFK in terms of the period andfluctuations of the oscillations and the associated thermody-namic cost of generating such dynamics, which can be quan-tified using TUR. Finally, we will underscore the optimalityof glycolytic oscillations emerging from octameric PFK bycomparing it with the dynamics anticipated for hypotheticalnon-octameric PFK model.
THEORETICAL BACKGROUNDThermodynamic uncertainty relation for oscillatory processes
Biological processes are inherently stochastic due to fluctu-ations in the cellular environment. Suppose that time trajecto-ries of a dynamical process generated from a constant thermo-dynamic drive are available. Then, in order to gain knowledgeabout the dynamics from a given trajectory with higher preci-sion, it is required to analyze a longer trajectory; however,generation of a longer time trajectory incurs larger thermo-dynamic cost. According to the thermodynamic uncertaintyrelation (TUR), there are trade-offs between the total entropyproduction ∆ S tot ( t ) from the process generated under a con-stant drive and the precision of a time integrated current-like observable X ( t ) to probe the dynamical process generated innonequilibrium. Furthermore, the product between ∆ S tot ( t ) and the squared relative uncertainty, Var [ X ( t )] / (cid:104) X ( t ) (cid:105) , de-fined as the uncertainty product Q , is bounded below by twice a r X i v : . [ phy s i c s . b i o - ph ] F e b the Boltzmann constant ( k B ): Q ≡ ∆ S tot ( t ) Var [ X ( t )] (cid:104) X ( t ) (cid:105) ≥ k B . (1)Because of the lower bound, there is a minimal thermody-namic (entropic) cost to generate a process with a certainprecision, or that the precision of the process is bounded bythe thermodynamic cost being expended. The relation ap-plies to current-like output observables, with odd parity un-der time reversal satisfying X ( t ) = − X ( − t ) [18], that arisefrom dynamic processes that can be modeled using continuoustime Markovian dynamics in discrete network or in continu-ous space under time-independent driving. There have beenfurther generalizations of the relation [17–24] including thethose under time-dependent driving [24–26]. Here, we con-fine ourselves to the original relation (Eq.1) since the oscil-latory dynamics considered in this study are generated undertime-independent driving. Throughout the paper, we set theBoltzmann constant k B to unity for simplicity of the expres-sions.There have been several applications of TUR to charac-terize biological processes [27–32]. For biological motorsthat move along filaments, the net displacement of a motor ∆ x ( t ) = x ( t ) − x (0) = − ∆ x ( − t ) = x (0) − x ( − t ) satis-fies the odd parity, with time translation symmetry, and maybe selected as a proper output observable to probe the func-tional dynamics of the motor. For dynamical processes dis-playing temporal oscillations, the time duration of the process, t , can be employed as a current-like observable to monitor theprogress of biochemical oscillations, i.e., X ( t ) = t . Then,the time duration, t , is decomposed as t = (cid:80) ni =1 T i where n is the number of oscillations for time t and T i is a stochasticvariable representing the period of i -th oscillation. Then, thefirst and second moments of the time duration are written as (cid:104) t (cid:105) = n (cid:88) i =1 (cid:104) T i (cid:105) = n (cid:104) T (cid:105) (2)and (cid:104) t (cid:105) = n (cid:88) i =1 (cid:104) T i (cid:105) + 2 (cid:88) i>j (cid:104) T i T j (cid:105) = n (cid:104) T (cid:105) + n ( n − (cid:104) T i T j (cid:105) = n (cid:104) T (cid:105) − n (cid:104) T (cid:105) + (cid:104) t (cid:105) (3)where (cid:104) T (cid:105) ≡ (1 /n ) (cid:80) ni =1 T i is the mean oscillatory period,and the independence of two oscillatory periods, (cid:104) T i T j (cid:105) = (cid:104) T i (cid:105)(cid:104) T j (cid:105) = (cid:104) T (cid:105) , are used in the last line of Eq.3. Thus, fromEqs. 2 and 3, one obtains (cid:104) t (cid:105)−(cid:104) t (cid:105) = n (cid:104) T (cid:105)− n (cid:104) T (cid:105) , and theuncertainty product Q for processes demonstrating temporal oscillations as follows [32–34]: Q = ∆ S tot ( t ) (cid:104) t (cid:105) − (cid:104) t (cid:105) (cid:104) t (cid:105) = (cid:20) ∆ S tot ( t ) n (cid:104) T (cid:105) (cid:21) Var ( T ) (cid:104) T (cid:105) = ˙ S tot Var ( T ) (cid:104) T (cid:105) = ∆ S cyc Var ( T ) (cid:104) T (cid:105) ≥ k B , (4)where Var ( T ) = (cid:104) ( δT ) (cid:105) = (cid:104) T (cid:105) − (cid:104) T (cid:105) . ˙ S tot =∆ S tot ( t ) /n (cid:104) T (cid:105) = ∆ S tot ( t ) /t is the rate of entropy produc-tion, and the entropy production per cycle can be defined as ∆ S cyc ≡ ∆ S tot ( t ) /n = ˙ S tot × (cid:104) T (cid:105) . The inequality in the lastline of Eq.4 is identical to the expression N ≡ (cid:104) T (cid:105) /D ≤ ∆ S cyc / ( k B = 1 ) of Marsland et al. ’s [32, 35] with D ≡ Var ( T ) = (cid:104) ( δT ) (cid:105) . If the uncertainty product Q evaluated fora process displaying temporal oscillations is close to the the-oretical minimum k B , one could presume that the process isclose to its optimal condition where the thermodynamic costof generating an oscillatory dynamics with a certain temporalprecision is minimal.For a given oscillatory process arising from a set of ki-netic equations, it is straightforward to evaluate the mean( (cid:104) T (cid:105) ) and variance of oscillation period (Var ( T ) ) in Eq.4 fromthe time trajectories of dynamical variables, x ( t ) , that satisfy x ( t ) (cid:39) x ( t + T ) . To obtain the total entropy production ratefrom x ( t ) , one can utilize the evolution equation of proba-bility density, P ( x , t ) , namely the Fokker-Planck (FP) equa-tion. The FP equation is obtained from a corresponding set ofchemical master equations (CME) at a finite volume ( Ω ) (seeSI): ∂P ( x , t ) ∂t = − (cid:126) ∇ · (cid:104) H · (cid:126)F ( x ) − D ( x ) · (cid:126) ∇ (cid:105) P ( x , t )= − (cid:126) ∇ · (cid:126)J ( x , t ) , (5)where D ( x )(= k B T H ( x )) is the Ω -dependent diffusion ten-sor with H ( x ) being the motility tensor, and (cid:126)J ( x , t ) = H ( x ) · (cid:126)F ( x ) P ( x , t ) − D ( x ) · (cid:126) ∇ P ( x , t ) (6)is the probability current, (cid:126)F ( x , t ) is the driving force vec-tor. The total entropy production rates, contributed by boththe system and reservoir, are obtained by averaging the corre-sponding trajectory-based entropy production rates over theprobability density P ( x , t ) , ˙ S tot ( t ) = (cid:104) ˙ s tot ( t ) (cid:105) = (cid:90) ˙ s tot ( t ) P ( x ( t ) , t ) d x , (7)where ˙ s tot ( t ) = ˙ s sys ( t ) + ˙ s res ( t ) with s sys ( t ) = − log P ( x ( t ) , t ) and ˙ s res ( t )(= (cid:126)F ( x , t ) · ˙ x ) /T r where T r isthe temperature of the reservoir [36–38]. Together with theexpression of probability current (Eq.6), Eq.7 yields [36] ˙ S tot ( t ) = (cid:90) (cid:126)J (cid:124) ( x , t ) · D − ( x ) · (cid:126)J ( x , t ) P ( x , t ) d x . (8)The entropy production rate at steady state, obtained from J ss ( x ) and P ss ( x ) , allows us to evaluate the uncertainty prod-uct Q (Eq.4). Brusselator
The Brusselator, a model for autocatalytic reactions thatdisplay sustained oscillatory dynamics in certain range of pa-rameters, offers all the ingredients of oscillatory dynamics re-quired to learn the dissipation and precision of glycolytic os-cillations to be studied.The model forms an open system, consisting of two chem-ical compounds X and Y reacting each other and being de-pleted while other source compounds A and B are constantlysupplied to the system at fixed concentrations. The reactionscheme for Brusselator is A k −→ X,B k −→ Y, X + Y k −→ X,X k −→ φ. (9)Properly non-dimensionalized (see SI), the rate equations atthe limit of an infinite volume (Ω → ∞ ) is given as dxdt = f ( x, y ) = a − x + x y,dydt = g ( x, y ) = b − x y, (10)where x , y , a , and b are the non-dimensionalized concen-trations of chemical species X , Y , A and B , respectively.The equations for the concentrations of x and y are nonlin-ear; thus when x and y are expanded around a fixed point ( x ∗ , y ∗ ) that satisfies f ( x ∗ , y ∗ ) = 0 and g ( x ∗ , y ∗ ) = 0 suchthat x = x ∗ + δx and y = y ∗ + δy , it yields a set of linearizedequations δ ˙ x = J ( x ∗ , y ∗ ) · δ x . (11)where δ x ≡ ( δx, δy ) (cid:124) and J ( x ∗ , y ∗ ) = (cid:20) f x f y g x g y (cid:21) is the Jacobian matrix evaluated at ( x ∗ , y ∗ ) = ( a + b, b/ ( a + b ) ) , where f x ≡ ∂ x f ( x, y ) | ( x,y )=( x ∗ ,y ∗ ) . Since the fluctu-ation δ x is expected to change with time as δ x ∼ e λt , thestability of the fixed point is determined by the real partsof eigenvalues ( λ ) obtained from the characteristic equation det [ λ I − J ( x ∗ , y ∗ )] = 0 where I is the identity matrix. Theeigenvalues are: λ = 12 ( τ ± (cid:112) τ − (12) with τ = − ( a + b ) − ( a − b ) / ( a + b ) and ∆ = ( a + b ) .Since ∆ > for any a and b , the sign of Re ( λ ) is deter-mined entirely by the sign of τ . The condition of τ = 0 ,leading to b + 3 ab + (3 a − b + a + a = 0 , determinesthe phase boundary (Fig.1A). The set of parameters ( a, b ) be-longing to the blue region of phase diagram (Fig.1A) leads to τ < , then the fixed points are stable, and the time trajectoryof ( x ( t ) , y ( t )) converges to ( x ∗ , y ∗ ) (see the lower rightmostpanel of Fig.1A). On the other hand, the region colored in red( τ > ) yields unstable fixed points that produce limit cycles(see Fig.S1A and B for the vector fields ( ˙ x ( t ) , ˙ y ( t )) depictedfor ( a, b ) leading to unstable and stable fixed points). For thecase of 2D phase plane, existence of limit cycles is alwaysguaranteed by the Poincar´e-Bendixson theorem [39].For a system with a finite volume, Ω , the FP equation de-rived from CME [40] for the Brusselator is obtained with thedriving force vector (cid:126)F ( x, y ) and diffusion tensor D ( x, y ) . (cid:126)F ( x, y ) = (cid:20) a − x + x yb − x y (cid:21) + 12Ω (cid:20) − / − xy + x / xy − x / (cid:21) , D ( x, y ) = 12Ω (cid:20) a + x + x y − x y − x y b + x y (cid:21) . (13)The ( a, b ) -dependent steady state probability distribution P ss ( x, y ) , obtained from an ensemble of trajectories generatedusing the Gillespie simulations, allows us to calculate the (cid:126)J ss (Eq.5), and hence ˙ S tot at the steady state using Eq.8.In the parameter range of ( a, b ) yielding unstable fixedpoints ( τ > ), we have computed ( a, b ) -dependent 2D di-agrams of various quantities at Ω = 1600 : ˙ S tot , (cid:104) T (cid:105) , (cid:104) δT (cid:105) ,amplitude of oscillations ( A ), and the integral current J cycle = (cid:72) (cid:126)J ss · d(cid:126)l (cid:46) (cid:72) d(cid:126)l (Fig 1). The entropy production rate ˙ S tot showoverall positive correlation with (cid:104) T (cid:105) , (cid:104) δT (cid:105) , and A , but notwith J cycle . The larger J cycle signifies faster oscillations. Itis noteworthy that the correlation or anti-correlation betweenthe quantities calculated here is relatively clear for the Brus-selator, but same is not necessarily true for glycolytic oscilla-tor (compare Fig.S2A and B). The (cid:104) δT (cid:105) , displaying a non-monotonic variation, is minimized at a basin of parameterspace around ( a, b ) ≈ (0 . , . . The product of ˙ S tot and (cid:104) δT (cid:105) / (cid:104) T (cid:105) gives rise to the 2D diagram of Q ( a, b ) (Fig.1),indicating that Q is minimized in the vicinity of the phaseboundary ( a, b ) = (0 . , . to Q ≈ .Importantly, Q is independent of the system size Ω . Theentropy production is an extensive quantity that linearly in-crease with Ω for a given time interval t . Thus, the entropyproduction rate scales with the volume as ˙ S tot ∼ Ω . Next, thefluctuation of the oscillatory period (cid:104) ( δT ) (cid:105) is proportional tothe magnitude of the Ω -dependent diffusion tensor D , suchthat (cid:104) ( δT ) (cid:105) ∼ D ∼ Ω − as defined in Eq.13, whereas (cid:104) T (cid:105) is decided independently from Ω . Taken together, the uncer-tainty product Q is a quantity independent of Ω , which canalso be confirmed using a host of simulations with varying Ω (see Fig.S3). FIG. 1. A . Dynamic phase diagram of the Brusselator. The steady state oscillation of limit cycle occurs in the pink region, left side of thephase 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 phaseplanes with the vector fields are depicted in Fig.S1. B . The 2D diagram of Q ( a, b ) corresponding to the phase region of limit cycles. C . The2D diagrams for other quantities, ˙ S tot , (cid:104) T (cid:105) , (cid:104) ( δT ) (cid:105) , J cycle , and A , as a function of a and b . RESULTS
The allosteric regulation of PFK1 and its substrate andproduct concentration-dependent enzymatic activity can beformulated using the strategy of Mono-Wyman-Changeauxmodel [48, 49]. The enzyme PFK1, a oligomeric complexconsisting of n catalytic and n regulatory sites to whichthe substrate and product, respectively, can bind, is equili-brated between tense (inactive, T ) and relaxed (active, R )states, which differ in terms of their conformations and bind-ing affinities to substrate and product (see Fig.2). T statecan only accommodate substrate molecules, and the subscript i (= 0 , , . . . n ) in T i denotes the number of substrates boundto the binding sites. On the other hand, R state can accommo-date both substrate and product molecules; the two subscripts i (= 0 , , . . . , and j (= 0 , , . . . , n ) of R ij denote the num-ber of substrates and products bound to the catalytic and reg-ulatory sites, respectively. In the absence of substrate, twoapo states of PFK1, R ( ≡ R ) and T states, are chemicallyequilibrated with the ratio, L = [ T ] / [ R ] called allostericconstant, which determines the degree of cooperativity of theenzyme. The glucose converting into substate S (fructose-6-phosphate, F6P) via multiple steps is injected at a constantrate ν 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 . Theincrease of F6P as a result of the catalytic processes of glycol-ysis in turn increases the amount of FBP, which positively reg-ulates the catalytic activity of PFK1 by promoting the T -to- R transitions. The nonlinear response of FBP-binding inducedautocatalytic activation of PFK1 generates the glycolytic os-cillation [50, 51].We assume that the binding and unbinding of substrate ( S )and product ( P ) to and from each binding site of R state occurwith the rates k SR, on , k SR, off , and k PR, on , k PR, off , and that only the S can bind/unbind to T state with k ST, on , k ST, off . Then, theconcentration of each state of PFK1 is obtained as follows byassuming a quasi-steady-state approximation [44]. [ T i ] = (cid:18) ni (cid:19) (cid:18) [ S ] K ST (cid:19) i [ T ] = L (cid:18) ni (cid:19) (cid:18) c [ S ] K SR (cid:19) i [ R ] , [ R ij ] = (cid:18) ni (cid:19) [ α ([ S ])] i (cid:18) nj (cid:19) [ γ ([ P ])] j [ R ] (14)where α ([ S ]) = [ S ]( K SR + k/k SR, on ) ,γ ([ P ]) = [ P ] K PR , FIG. 2. Allosteric model for PFK in glycolysis. The PFK1 is an oligomeric enzyme, which takes different oligomeric state depending on theorganism. Depicted is the octameric structure (dimer of tetramer) of PFK1 in R state (PDB code 4U1R [41]). 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 modelof glycolytic oscillation assumes an open thermodynamic system where glucose (or substrate) molecules are injected at the rate ν and theproduct ( P ) is drained at k s . Depicted is the octameric form of enzyme with n = 8 binding sites for substrate and product. At low substrateconcentration, the enzyme is in the tense ( T ) state. But, with increasing substrate concentration, transition to the R state occurs, increasing thecatalytic activity of enzyme. K SR (= k SR, off /k SR, on ) and K ST (= k ST, off /k ST, on ) are the bindingaffinities (dissociation constants) of the substrate to the cat-alytic site in the R and T states, respectively, whereas K PR is the binding affinity of the product to the regulatory site inthe R state. The parameter c = K SR /K ST denotes the ratio ofthe dissociation constants of the substrate from the catalyticsites in R and T states. Then, from the concentration of eachenzyme state [ T i ] and [ R ij ] (Eq.14), it is straightforward tocalculate a binding polynomial ¯ Y ( ≤ ¯ Y ≤ ), namely thefraction of catalytic sites in the R state bound by the substrate, ¯ Y = total substrates bound to R state n × total enzymes = (cid:80) ni =0 (cid:80) nj =0 i [ R ij ] nZ ([ S ] , [ P ])= α ([ S ]) q ([ S ] , [ P ]) nZ ([ S ] , [ P ]) (15) where q ([ S ] , [ P ]) ≡ n (1 + α ([ S ])) ( n − (1 + γ ([ P ])) n [ R ] , and the total concentration of enzyme Z ([ S ] , [ P ]) Z ([ S ] , [ P ]) = n (cid:88) i =0 [ T i ] + n (cid:88) i =0 n (cid:88) j =0 [ R ij ]= (cid:20) L (cid:18) c [ S ] K SR (cid:19) n + (1 + α ([ S ])) n (1 + γ ([ P ])) n (cid:21) [ R ] (16)Substrates supplied with constant rate ν bind to T or R state, and is catalyzed by the R state of PFK1 ( R ij ) with arate k to generate product molecules, whereas the products ei-ther bind to the regulatory site of R state or depleted from thesystem at a rate k s . The mass action laws of the substrate and FIG. 3. A . Phase diagram of glycolytic oscillations as a function of ν and k s with the allosteric constant ( L = 4 × ), catalytic rate ( k = 500 s − ) and c = 0 . . No catalysis occurs in the hashed region. The simulations were performed with ν = 0 . − . mM/s, demarcated withthe dark pink region, corresponding to the glucose uptake rate of yeast. B . Three representative trajectories showing oscillatory dynamics ofsubstrate and product concentrations. The trajectory marked with yellow star was generated at the parameter values relevant for octameric PFKof yeast giving rise to the oscillatory period of ∼
400 s [42–47] ( L = 4 × , k = 500 s − , ν = 0 . mM/s, k s = 0 . s − , c = 0 . ). C .2D diagram of uncertainty product, Q as a function of ν and k s . D . ˙ S tot , (cid:104) T (cid:105) , (cid:104) ( δT ) (cid:105) , J cycle , and A calculated as a function of ν and k s . product yield a set of coupled nonlinear equations [44]: d [ S ] dt = ν − k n (cid:88) i =0 n (cid:88) j =0 i [ R ij ] = ν − kα ([ S ]) q ([ S ] , [ P ]) ,d [ P ] dt = k n (cid:88) i =0 n (cid:88) j =0 i [ R ij ] − k s [ P ]= kα ([ S ]) q ([ S ] , [ P ]) − k s [ P ] . (17)Here, by specifically considering the yeast PFK which adoptsan octameric form ( n = 8 ), we use the allosteric constant L =4 × which is 3 orders of magnitude greater than the valueknown for the tetramer ( n = 4 ) [42]. Other parameters relatedwith substrate/product binding to each protomer are expectedto be identical between tetramer and octamer, thus we take k SR, on = 2000 /mM · s, K SR = 0 . mM and K PR = 0 . mM[42, 45–47]. The ratio of the substrate binding affinities to the R and T states is c = K R /K T = 0 . , so that the substratebinds preferentially to the R state. The ATP hydrolysis timedue to ATPase activity is typically (cid:38) O (1) msec [52], andhence we set k = 500 s − .Following the same procedure used in the analysis of Brus-selator, we analyze Eq.17 to calculate phase diagrams as afunction of ν and k s (Fig.3A) and generate dynamical trajec-tories of [ S ]( t ) and [ P ]( t ) (Fig.3B). For the values of ( ν, k s ) pertaining to the limit cycles, the substrate concentration dis-play saw-tooth like oscillatory pattern in time, and and theproduct concentration spikes when the substrate concentrationfalls, which generates a loop in the phase plane of ([ S ] , [ P ]) (Fig.3B). A trajectory with the mean oscillatory period of ∼
400 s ( ≈ − min), fluctuating between ∼ ∼ L = 4 × , ν = 0 . mM/s, k s = 0 . s − , c = 0 . ) (yellow stars in Figs. 3A,B, C). The period and amplitude of the oscillation comportwell with those observed in yeast and yeast extract [43, 53], inwhich PFK1 enzymes are oligomerized to an octameric form.Next, the FP equation for the glycolytic oscillations is ob-tained from the CME corresponding to Eq.17 with the follow-ing force vector and diffusion tensor: (cid:126)F = (cid:20) ν − kαqkαq − k s [ P ] (cid:21) + 12Ω (cid:20) − kαq − kα∂ [ S ] q + kα∂ [ P ] q − k s + kαq + kα∂ [ S ] q − kα∂ [ P ] q (cid:21) , D = 12Ω (cid:20) ν + kαq − kαq − kαq k s [ P ] + kαq (cid:21) , (18)where the concentration dependences of α = α ([ S ]) (Eq.15)and q = q ([ S ] , [ P ]) (Eq.16) are omitted for the simplicity ofthe expression. The FP equation with these (cid:126)F and D is used tocalculate J ss ([ S ] , [ P ]) , P ss ([ S ] , [ P ]) , and ˙ S tot based on Eq.8.Shown in Fig.3C, D are the 2D diagrams of ˙ S tot , (cid:104) T (cid:105) , (cid:104) δT (cid:105) , J cycle , A , and finally Q as a function of ν and k s , whichare the two experimentally controllable parameters. Overall,the correlations between these quantities are not so strongin comparison with those calculated for the Brusselator (seeFig.S2). It is fair to say that the correlation or a trend seenin Brusselator cannot be generalized to other biochemical os-cillators (see Fig.S2). At the parameter values, ν = 0 . mM/s − and k s = 0 . s − , yielding T = (cid:104) T (cid:105) ± (cid:104) δT (cid:105) / ≈ ± s, the uncertainty product is Q (cid:39) . Remark-ably, while Q (cid:39) is observed in the vicinity of the lowerphase boundary where the entropy production rate is minimalover the relevant phase space (see Fig.3C and the first panelof Fig.3D), (cid:104) T (cid:105) ≈ sec (the second panel of Fig.3D) isobtained only at ( ν, k s ) (cid:39) (0 . mM/s , . s − ) . DISCUSSIONS
Glycolytic oscillations are the temporal order that emergesunder certain special conditions in which parameters defin-ing a set of coupled nonlinear equations yield unstable fixedpoints. Notably, for octameric form of PFK oligomers todemonstrate oscillatory dynamics, the condition of c = K R /K T = 0 . , which renders the substrate binding tothe protomer in R state more preferable than to T state bya hundred fold, is essential. If c is increased to 0.1, thephase space corresponding to limit cycles significantly nar-rows down (Fig.4). For tetrameric form of PFK ( n = 4 ),which pertains to bacteria, the phase space region for limit cy-cles is much narrower even when c = 0 . and L is set to thevalue of octamer ( L = 4 × ) (Fig.4). Unless k s is tunedto a narrow interval of k s (cid:39) . − . s − , no oscillationis expected. The ( ν, k s ) phase diagrams of glycolysis withvarying c and n (Fig.4) rationalize why glycolytic oscillationswere only reported in eukaryotic PFK, where PFK exists inthe octameric form.The TUR, which specifies the physical lower bound to theuncertainty product, is used to assess how the period of tem-poral order emerging from the underlying dynamical processis balanced with the dissipation under the constraint of cost-precision trade-off. In the Brusselator, whose dynamical be-havior is defined only with two parameters ( a and b ), thereis a specific case that both precision of oscillatory periodand dissipation are simultaneously minimized to yield a rea-sonably small uncertainty product Q ≈ over the phasespace. In comparison, glycolytic oscillations are more com-plicated with many more parameters ( L , c , n , K SR , K PR , k , ν , k s ). To simplify the problem, we have reduced the un-knowns by assuming that some of the parameters have iden-tical values with those pertaining to the protomer. The val-ues of uncertainty product Q for the glycolytic oscillations attheir working condition producing the period of ∼ (5 – 10)min is Q (cid:39) . Remarkably, given the substrate injection rate ν = 0 . − . s − , Q (cid:39) is effectively the minimalvalue over the phase space involving the limit cycle (Fig.3C). Q (cid:39) is greater than those determined for the molecularmotors Q ≈ − [28, 54], and biological copy process byexonuclease-deficient T7 DNA polymerase Q ≈ [30], butsmaller than Q ≈ − for the translation process by E.coli ribosome [30, 31]. In comparison with the uncertaintyproduct determined for other biochemical cycles (
Q ≈ ) FIG. 4. The 2D phase diagram of glycolytic oscillations as a functionof the injection ( ν ) and degradation rates ( k s ) for varying c at n (=2 , , with L = 4 × and k = 500 s − . No catalysis occurs inthe hashed region. The working condition of glycolytic oscillations( ν = 0 . mM/s, k s = 0 . s − ) that produces the oscillationswith the mean period (cid:104) T (cid:105) = 400 s at c = 0 . and n = 8 (Fig.3) ismarked with the yellow star in each panel. [32], which severely underperform the TUR’s lower boundof , the value of the uncertainty product Q (cid:39) for theglycolytic oscillation arising from octameric PFK is signifi-cantly smaller, minimizing the entropy production rate overthe relevant phase, which indicates the cost-effectiveness ofthe molecular mechanism generating the oscillatory dynam-ics.Lastly, it is of particular note that the ∼ (5 – 10) min oscil-lation period is observed in cellular or physiological scales aswell, such as the blood glucose level, intracellular Ca con-centrations, and membrane action potentials, which are con-trolled by the pulsatile secretion of insulin with a period of ∼ (5 – 10) min [55–57], suggestive of a connection between thedynamics at the molecular and macroscopic scales [13], andtheir synchronization [58]. Remarkably, the loss of pulsatileinsulin release resulting from elevated glucose level [55, 58]is also consistent with our study that the working conditionof glycolytic oscillation is situated near the borderline of thephase boundary. The phase diagram depicted in Fig.3A pre-dicts that a moderate elevation of glucose injection rate be-yond ν ≈ . mM/s for a fixed k s = 0 . s − would abolishthe oscillations.We thank Prof. Junghyo Jo for illuminating discussions onglucose level oscillations. This study is supported by KIASIndividual Grants CG076501 (P.K.) and CG035003 (C.H.).We thank the Center for Advanced Computation in KIAS forproviding computing resources. ∗ [email protected][1] C. Bustamante, J. Liphardt, and F. Ritort, Physics Today , 43(2005).[2] M. J. Rust, S. S. Golden, and E. K. O’Shea, Science , 220(2011).[3] M. J. Rust, J. S. Markson, W. S. Lane, D. S. Fisher, and E. K.O’Shea, Science , 809 (2007).[4] B. Novak and J. Tyson, J. Cell Sci. , 1153 (1993).[5] J. Tyson and B. Novak, J. Theor. Biol. , 249 (2001).[6] J. E. Ferrell, T. Y.-C. Tsai, and Q. Yang, Cell , 874 (2011).[7] I. Prigogine, Science , 777 (1978).[8] A. Goldbeter, Phil. Trans. Roy. Soc. A: Math. Phys. Eng. Sci. , 20170376 (2018).[9] L. Duysens and J. Amesz, Biochimi. Biophys. Acta , 19(1957).[10] A. Boiteux, A. Goldbeter, and B. Hess, Proc. Natl. Acad. Sci.U. S. A. , 3829 (1975).[11] A. Goldbeter and S. R. Caplan, Annu. Rev. Biophys. Bioeng. ,449 (1976).[12] K. Tornheim, Diabetes , 1375 (1997).[13] R. Bertram, L. Satin, M. Zhang, P. Smolen, and A. Sherman,Biophys. J. , 3074 (2004).[14] B. Hess and A. Boiteux, Annu. Rev. Biochem. , 237 (1971).[15] A. C. Barato and U. Seifert, Phys. Rev. Lett. , 158101(2015).[16] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England,Phys. Rev. Lett. , 120601 (2016).[17] J. M. Horowitz and T. R. Gingrich, Nat. Phys. , 15 (2019).[18] Y. Hasegawa and T. Van Vu, Phys. Rev. Lett. , 110602(2019).[19] K. Proesmans and C. Van den Broeck, EPL , 20001 (2017).[20] J. S. Lee, J.-M. Park, and H. Park, Phys. Rev. E , 062132(2019).[21] B. K. Agarwalla and D. Segal, Phys. Rev. B , 155438 (2018).[22] S. Lee, C. Hyeon, and J. Jo, Phys. Rev. E , 032119 (2018).[23] P. P. Potts and P. Samuelsson, Phys. Rev. E , 052137 (2019).[24] T. Koyuk, U. Seifert, and P. Pietzonka, J. Phys. A: Math. Theor. , 02LT02 (2018).[25] T. Koyuk and U. Seifert, Phys. Rev. Lett. , 230601 (2019).[26] T. Koyuk and U. Seifert, Phys. Rev. Lett. , 260604 (2020).[27] P. Pietzonka, A. C. Barato, and U. Seifert, J. Stat. Mech. TheoryExp. , 124004 (2016).[28] W. Hwang and C. Hyeon, J. Phys. Chem. Lett. , 513 (2018).[29] M. Uhl and U. Seifert, Phys. Rev. E. , 022402 (2018).[30] Y. Song and C. Hyeon, J. Phys. Chem. Lett. , 3136 (2020).[31] W. D. Pi˜neros and T. Tlusty, Phys. Rev. E , 022415 (2020).[32] R. Marsland III, W. Cui, and J. M. Horowitz, J. Roy. Soc. In-terface , 20190098 (2019).[33] Y. Cao, H. Wang, Q. Ouyang, and Y. Tu, Nature Physics (2015), 10.1038/nphys3412.[34] L. G. Morelli and F. J¨ulicher, Phys. Rev. Lett. , 228101(2007).[35] R. Marsland III and J. England, Rep. Prog. Phys. , 016601(2017).[36] U. Seifert, Phys. Rev. Lett. , 040602 (2005).[37] H. Ge and H. Qian, Phys. Rev. E , 051133 (2010).[38] T. Tom´e and M. J. de Oliveira, Phys. Rev. E , 021120 (2010).[39] S. H. Strogatz, Nonlinear dynamics and chaos with student so-lutions manual: With applications to physics, biology, chem-istry, and engineering (CRC press, 2018).[40] H. Qian, S. Saffarian, and E. L. Elson, Proc. Natl. Acad. Sci. U. S. A. , 10376 (2002).[41] M. Kloos, A. Br¨user, J. Kirchberger, T. Sch¨oneberg, andN. Str¨ater, Biochem. J. , 421 (2015).[42] D. Blangy, H. Buc, and J. Monod, J. Mol. Biol. , 13 (1968).[43] B. Hess, A. Boiteux, and J. Kr¨uger, Adv. Enzyme Reg. , 149(1969).[44] A. Goldbeter and R. Lefever, Biophys. J. , 1302 (1972).[45] S. M. M. McCoy, M. A. and D. F. Wyss, J. Am. Chem. Soc. , 149 (2005).[46] J. W. Cardon and P. D. Boyer, Eur. J. Biochem. , 443 (1978).[47] M. Moyer, S. Gilbert, and K. Johnson, Biochemistry ,800—813 (1998).[48] J. Mono, J. Wyman, and J. P. Changeux, J. Mol. Biol. , 88(1965).[49] D. Thirumalai, C. Hyeon, P. I. Zhuravlev, and G. H. Lorimer,Chem. Rev. , 6788 (2019).[50] K. Tornheim, J. Biol. Chem. , 2619 (1988).[51] G. C. Yaney, V. Schultz, B. A. Cunningham, G. A. Dunaway,B. E. Corkey, and K. Tornheim, Diabetes , 1285 (1995).[52] S. P. Gilbert and K. A. Johnson, Biochemistry , 1951 (1994).[53] V. Hess and A. Szabo, J. Chem. Edu. , 289 (1979).[54] M. L. Mugnai, C. Hyeon, M. Hinczewski, and D. Thirumalai,Rev. Mod. Phys. , 025001 (2020).[55] J. P. McKenna, R. Dhumpa, N. Mukhitov, M. G. Roper, andR. Bertram, PLoS Comp. Biol. , e1005143 (2016).[56] D. A. Lang, D. R. Matthews, J. Peto, and R. C. Turner, NewEngland Journal of Medicine , 1023 (1979).[57] P. O. Westermark and A. Lansner, Biophys. J. , 126 (2003).[58] B. Lee, T. Song, K. Lee, J. Kim, S. Han, P.-O. Berggren, S. H.Ryu, and J. Jo, PLoS one , e0172901 (2017). SUPPORTING INFORMATIONNon-dimensionalization of rate equations
The stochastic version of the Brusselator is written as dXdt = k A − k X + k X ( X − Y Ω ,dYdt = k B − k X ( X − Y Ω (S1)where Ω is the volume of the system. Using the followingtransformations of the variables and parameters, (cid:18) k k (cid:19) / X Ω −→ x, (cid:18) k k (cid:19) / Y Ω −→ y,k t −→ t,k k (cid:18) k k (cid:19) / A Ω −→ a,k k (cid:18) k k (cid:19) / B Ω −→ b, (S2)one can write down a non-dimensionalized version of the rateequations in Eq.10 at the limit of Ω → ∞ . Fokker Planck equation from Chemical Master Equations
For a chemical species (X) involved in a reaction: sX + · · · k −→ s (cid:48) X + · · · (S3)the time evolution equation of X at the deterministic limit canbe written as d [ X ] dt = Sk [ X ] s (S4)where [ X ] is the concentration of X , S = s (cid:48) − s is thestoichiometric coefficient. The chemical state of the systemat any time is fully determined by the state vector X =( X , . . . , X N ) where X i is the number of chemical species X i in a compartment of a finite volume Ω . Then, probabilitydistribution for the system to be in state X at time t is P ( X , t + dt ) = P ( X , t )+ dt R (cid:88) r =1 [ f r ( X − S r ) P ( X − S r , t ) − f r ( X ) P ( X , t )] . (S5) where f r ( X ) = k r Ω (cid:81) Ni =1 X i !( X i − s ir )!Ω sir is the probability forthe reaction r to occur. At the limit of dt → , Eq.S5 is castinto the Chemical Master Equation (CME) ∂ t P ( X , t ) = R (cid:88) r =1 { f r ( X − S r ) P ( X − S r , t ) − f r ( X ) P ( X , t ) } . (S6)Since analytic solutions of CME is known only for limitedcases, CME is typically approximated to Fokker Planck equa-tion through a Taylor expansion of the relevant terms to thesecond order, f r ( X − S r ) P ( X − S r , t ) − f r ( X ) P ( X , t ) ≈ − N (cid:88) i =1 S ir ∂ X i [ f r ( X ) P ( X, t )]+ N (cid:88) i,j =1 S ir S jr ∂ X i ∂ X j [ f r ( X ) P ( X, t )] , (S7)leading to ∂ t P ( X , t ) = − N (cid:88) i =1 ∂ X i [ A i ( X ) P ( X, t )]+ 12 N (cid:88) i,j =1 ∂ X i ∂ X j [ B ij ( X ) P ( X, t )] (S8)where the drift vector A and diffusion matrix B are given by A i ( X ) = R (cid:88) r =1 S ir f r ( X ) ,B ij ( X ) = R (cid:88) r =1 S ir S jr f r ( X ) . (S9)In order to handle the time evolution of different chemi-cal species inside the compartment of volume Ω in termsof their concentrations, we define x i ≡ X i / Ω ≡ [ X i ] .With this definition, the probability density of a set of con-centrations can be converted to that of molecular counts as P ( x , t ) ≡ Ω N P (Ω x , t ) . Finally, the Fokker-Planck equationfor the stochastic time evolution of concentrations of chemicalspecies is obtained as ∂ t P ( x , t ) = − N (cid:88) i =1 ∂ x i [ A i ( x ) P ( x, t )]+ 12Ω N (cid:88) i,j =1 ∂ x i ∂ x j [ B ij ( x ) P ( x, t )] . (S10)0 FIG. S1. Vector field ( dx/dt, dy/dt ) of Brusselator around the fixedpoint depicted on the phase plane. (A) Unstable fixed point with a = 0 . and b = 0 . . (B) Stable fixed point with a = 0 . and b = 0 . . y = 1 /x − a/x and y = b/x are drawn in pink and paleblue, respectively. · S tot · S tot ⟨ T ⟩⟨ δT ⟩ J cycle A A B