{"paper_id":"1f3bde0a-b61a-4012-8248-bca71e0e8618","body_text":"1 \nAnalysis of Metabolic Stability under Environmental Perturbations Using a Kinetic Model 2 \n 3 \nAtsuki Hishida1, Y usuke Himeoka2*, Chikara Furusawa2 3*  4 \n 5 \n1. Graduate School of Science, Kyoto University, Kyoto, Japan 6 \n2. Universal Biology Institute, The University of Tokyo, Tokyo, Japan 7 \n3. RIKEN Center for Biosystems Dynamics Research, Kobe, Japan 8 \n 9 \n*Corresponding author 10 \nE-mail: yhimeoka@ubi.s.u-tokyo.ac.jp (YH), furusawa@ubi.s.u-tokyo.ac.jp (CF) 11 \n 12 \nAbstract 13 \nCells robustly maintain metabolic function despite environmental fluctuations affecting reaction 14 \nkinetics, yet kinetic models often exhibit fragility. We investigated this discrepancy by examining the 15 \neffects of temperature changes on a kinetic model of Escherichia coli central metabolism. A gradual 16 \ntemperature decrease caused an abrupt transition to a state with sharply reduced ATP production 17 \nefficiency, triggered by an elevated ATP/ADP ratio creating a glycolysis bottleneck. Introducing a 18 \nrapid ATP-ADP exchange reaction prevented this ATP/ADP ratio surge and sustained high ATP 19 \nproduction efficiency across temperatures. Furthermore, altering enzyme abundances similarly 20 \nmaintained high A TP production, with predicted enzyme regulation aligning with experimental 21 \nobservations in E. coli at low temperatures. Our findings indicate that balancing key cofactors, such 22 \nas the ATP/ADP ratio, is crucial for preserving metabolic stability under environmental perturbations. 23 \n 24 \n  25 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nIntroduction 26 \nCellular metabolism comprises thousands of reactions that together generate energy, build 27 \nbiomolecules, and recycle cellular components. As external conditions fluctuate, cells adjust their 28 \nmetabolic state to maintain function. For example, bacteria switch between fermentative and 29 \nrespiratory regimes, reroute carbon through the pentose phosphate pathway (PPP) under oxidative 30 \nstress, and upregulate anaplerotic reactions when nutrients are scarce (Christodoulou et al., 2018; 31 \nSawers, 2025; Sudarsan et al., 2025). This plasticity arises from regulatory mechanisms such as 32 \nenzyme reallocation and allosteric modulation, enabling cells to occupy distinct, environment-specific 33 \nmetabolic states while sustaining growth (Gruber et al., 2025; Miyakoshi, 2024). 34 \nMathematical models of metabolic dynamics have been widely used to probe this plasticity 35 \n(El-Samad et al., 2005; Stephanopoulos et al., 1998; Strutz et al., 2019). Such models allow analysis 36 \nof how steady states shift as parameters vary. A recurring finding is that steady-state stability is often 37 \nfragile to parameter changes. For instance, a kinetic model of E. coli core metabolism was stable for 38 \nonly a small fraction of flux and concentration sets consistent with its kinetic parameters (Chakrabarti 39 \net al., 2013). Similarly, Steuer et al. (2006) sampled thousands of parameter sets for a kinetic model 40 \nof yeast glycolysis and found that over 80% of resulting steady states were unstable. These narrow 41 \nstability ranges suggest that environmental perturbations can readily destabilize cellular metabolism. 42 \nYet, despite the fragility predicted by models, real cells maintain stable metabolic states as 43 \ntemperature and nutrient availability change (Phadtare, 2004; Phadtare & Inouye, 2008; Riccardi et 44 \nal., 2023; Zhang & Gross, 2021). This discrepancy implies the existence of mechanisms that stabilize 45 \nmetabolism. Indeed, cells adjust enzyme levels and pathway activities in response to environmental 46 \ncues (Jozefczuk et al., 2010; Schmidt et al., 2016; Ye et al., 2012; Zhao et al., 2024). However, how 47 \nsuch regulation preserves dynamic stability remains unclear. 48 \nAmong environmental perturbations, temperature shifts directly and broadly affect 49 \nmetabolism (Noll et al., 2020; Wendering & Nikoloski, 2023). Temperature modulates the rate 50 \nconstants of all intracellular reactions, with magnitudes set by each reaction’s activation free energy 51 \n(Evans & Polanyi, 1935; Eyring, 1935). Because forward and reverse reactions need not respond 52 \nidentically, temperature changes can alter both rates and reversibility, reshaping flux throughout the 53 \nnetwork. Such nonuniform effects may destabilize steady states. Clarifying how cells prevent this 54 \noutcome would reveal general principles of metabolic robustness. 55 \nHere, we investigate mechanisms that keep complex metabolic systems stable under 56 \ntemperature-induced parameter perturbations. Using a kinetic model of E. coli central metabolism 57 \nconstructed from experimental data, we varied model parameters continuously with temperature and 58 \nobserved a bifurcation: the system transitioned discontinuously to a new steady state. The post-59 \nbifurcation state exhibited an elevated A TP/ADP ratio and reduced A TP production efficiency. 60 \nImposing operations that maintain the ATP/ADP ratio during the temperature shift eliminated the 61 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nbifurcation and preserved ATP production efficiency even at low temperatures. Moreover, adjusting 62 \nenzyme abundances to increase A TP-utilizing reaction rates stabilized the ATP/ADP ratio. These 63 \nresults indicate that maintaining the balance of key cofactors is critical for stabilizing metabolic 64 \nsystems under parameter perturbations caused by environmental change. 65 \n 66 \nResults 67 \nTemperature decrease induces a bifurcation 68 \nTo examine how a low-temperature shift affects the metabolic state of E. coli central metabolism, we 69 \nsimulated metabolic dynamics using rate constants adjusted for different temperatures. Figure 1 shows 70 \nthe reaction network of the kinetic model used in this study (see Materials and Methods for details). 71 \nRate constants at each temperature were calculated using the Eyring–Polanyi equation. We focused on 72 \ntemperature downshifts because heating can cause irreversible enzyme inactivation, which is not 73 \ncaptured by the standard Eyring–Polanyi formulation (Arcus & Mulholland, 2020; Noll et al., 2020). 74 \nActivation free energies (∆𝐺‡) for individual reactions were obtained from parameter values fitted to 75 \nthe in vivo metabolic state at 37 °C reported by Khodayari et al. (2014) (see Materials and Methods 76 \nfor the details). Without a temperature shift, the system converges to a stable steady state in which 77 \nglucose serves as the sole carbon source and ATP is produced via glycolysis, the TCA cycle, and the 78 \nelectron transport chain (ETC). We refer to this as the “original steady state” throughout the paper. 79 \nWhen the temperature is lowered from 37 °C to 0 °C, each rate constant decreases according 80 \nto its activation free energy. As a result, the distribution of rate constants shifts, and the median value 81 \ndecreases by about 10-fold (Figure 2A). We simulated system dynamics using the rate constants at 82 \n0 °C, starting from the original steady state. Numerical simulations showed that metabolite 83 \nconcentrations converged to a new stable steady state (Figure 2B). Compared with the original state, 84 \nmetabolite concentrations in this new steady state changed by factors ranging from 10ିଽ  to 10ସ —85 \nfar greater than the fold change observed in the median rate constant. 86 \nTo determine the temperature at which this substantial shift occurs, we tracked the steady 87 \nstate as the temperature was gradually decreased from 37 °C to 0 °C. At each step, system dynamics 88 \nwere integrated up to 𝑡 = 10ଵ଴  𝑠 to ensure convergence, and the resulting concentrations were used 89 \nas initial conditions for the next step, with temperature reduced by 0.1 °C. At 14.2 °C, the steady state 90 \nshifted discontinuously, with abrupt changes in the concentrations of all metabolites (Figure 2C shows 91 \none example; Supplementary Figure S1 presents all metabolite concentrations at the converged state 92 \nacross temperatures). Reaction rates also changed discontinuously, and relative changes compared 93 \nwith the original steady state are shown in Figure 2D.  94 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\n 95 \n  96 \nFigure 1. Reaction network of the kinetic model used in our simulations. Although each reaction was \ndecomposed into elementary steps for the simulations, the network is shown here with the pre -\ndecomposition reactions for clarity. \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nThis discontinuous shift could arise from either the disappearance or destabilization of the 97 \noriginal steady state, or from a switch between coexisting attractors. To discriminate between these 98 \npossibilities, we enumerated multiple stable steady states across temperatures including the transition 99 \npoint. We generated 10,000 random initial conditions by scaling each chemical concentration of the 100 \noriginal steady state by 10௎ (ିଷ .଴ ,ଷ .଴ ). Here 𝑈(𝑎, 𝑏) denotes a uniform distribution ranging from 𝑎 101 \nto 𝑏. We sampled initial conditions uniformly in the logarithmic scale over the interval (10ିଷ , 10ଷ ). 102 \nSee Materials and Methods for details of sampling procedure. Starting from each of these points, we 103 \nnumerically integrated the system dynamics to 𝑡 = 10ଵ଴  𝑠 and checked convergence for all cases.  104 \nThis procedure revealed three coexisting stable steady states above 14.2 °C and two below 105 \nit. These results indicate that the discontinuous shift arises from a bifurcation, in which the 106 \nphysiologically relevant steady state above 14.2 °C becomes destabilized as temperature decreases. 107 \nBifurcation diagrams for all small molecules are provided in Supplementary Figure S1. Notably, the 108 \ndisappearing steady state corresponds to the physiological E. coli  state, characterized by aerobic 109 \nrespiration and A TP production from glucose. 110 \nExamination of the steady state after the bifurcation revealed substantial changes in 111 \nmetabolic fluxes and ATP production efficiency. The discontinuous shift caused glycolytic and TCA 112 \ncycle fluxes to decrease by factors of 100–1,000, while flux through the pentose phosphate pathway 113 \n(PPP) increased (Figure 2D). As a result, ATP yield—defined as the sum of ATP-producing fluxes 114 \ndivided by the glucose uptake rate—declined following the temperature-induced bifurcation (Figure 115 \n2E). 116 \nTo investigate the cause of the decreased ATP yield, we examined metabolite concentrations 117 \nalong glycolysis and the TCA cycle. While most metabolite levels decreased after the temperature-118 \ninduced bifurcation, the concentration of D-glycerate 1,3-diphosphate (13dpg) increased by two orders 119 \nof magnitude compared with its pre-bifurcation level (Figure 2C; Supplementary Figure S1). Despite 120 \nthis elevation, the reaction degrading 13dpg—catalyzed by phosphoglycerate kinase (PGK), which 121 \nconverts 13dpg and ADP to 3-phosphoglycerate and ATP in glycolysis (Figure 1)—showed a reduced 122 \nrate. This reduction was caused by a ~40-fold increase in the ATP/ADP ratio, driven by elevated ATP 123 \nand decreased ADP concentrations (Figure 2F; Supplementary Figure S1). Because the PGK reaction 124 \nlies near the center of glycolysis, its slowdown created a bottleneck that diminished flux through both 125 \nglycolysis and the downstream TCA cycle.  126 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\n 127 \n  128 \nFigure 2. (A) Distribution of reaction rate constants in the system. (B) Time courses of relative \nmetabolite concentrations with all rate constants fixed at their values for 0ௗ°C and initial concentrations \nset to the original steady state; each line represents an individual metabolite. (C) Steady- state \nconcentration of D-glycerate 1,3- bisphosphate at each temperature. (D) Relative reaction rates at \nsteady state for each temperature, normalized to the rates at the original steady state. (E, F) Steady-\nstate ATP yield (E) and A TP/ADP ratio (F) at each temperature. \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nFixing the ATP/ADP ratio prevents the temperature-induced bifurcation 129 \nThe above results suggested that a sudden rise in the ATP/ADP ratio can trigger the temperature-130 \ninduced bifurcation. To investigate this relationship, we maintained the A TP/ADP ratio nearly constant 131 \nby introducing a rapid, reversible ATP–ADP exchange reaction into the system. Because the ATP/ADP 132 \nratio at the steady state before the bifurcation was approximately 1:1, the forward and reverse rate 133 \nconstants were set to 10⁶, far exceeding the median values of other reactions in the model. A detailed 134 \nprocedure for fixing the ATP/ADP ratio during temperature decrease is provided in Materials and 135 \nMethods. 136 \nThe modified model with the ATP–ADP exchange reaction maintained a single, stable 137 \nsteady state under the 37 ௗ°C parameter conditions. The resulting concentrations of 13dpg, the 138 \nATP/ADP ratio, and the ATP yield were similar to those of the original steady state (Figure ௗ3A–C). 139 \nWith the exchange reaction, lowering the temperature from 37ௗ°C to 0ௗ°C no longer caused a sharp rise 140 \nin 13dpg concentration or discontinuous changes in pathway fluxes, as shown by the relative changes 141 \nin reaction rates (Figureௗ3D). Consistently, the A TP/ADP ratio remained near unity, and the ATP yield 142 \nshowed no abrupt decline across the temperature range (Figureௗ3B,ௗC). To test for multiple attractors 143 \nat each temperature, we generated 10,000 random initial points by scaling each metabolite 144 \nconcentration of the original steady state by 10௎ (ିଷ .଴ ,ଷ .଴ )  and simulated the dynamics until 145 \nconvergence. All trajectories converged to a single steady state at every temperature. These results 146 \nindicate that, if the ATP/ADP ratio is fixed at approximately 1.0, the system does not undergo 147 \nbifurcation, even though individual reaction rate constants still vary with temperature. 148 \n  149 \nMaintaining the ATP yield by adjusting enzyme abundances 150 \nThe analysis in the previous section demonstrated that the ATP yield can be maintained at a high level 151 \nby fixing the ATP/ADP ratio using a rapid, reversible ATP–ADP exchange reaction. However, cells do 152 \nnot naturally possess such a reaction, and precise control of the ATP/ADP ratio in this manner is 153 \nbiologically implausible. Instead, cells adjust intracellular enzyme abundances to adapt to diverse 154 \nenvironmental perturbations, including temperature fluctuations, osmotic stress, oxidative conditions, 155 \nand nutrient limitations (Jozefczuk et al., 2010; Schmidt et al., 2016; Ye et al., 2012; Zhao et al., 2024). 156 \nWe therefore investigated whether ATP yield could be maintained under decreasing temperatures by 157 \nmodulating enzyme abundances. To this end, we sampled steady states with varied enzyme levels 158 \nunder low-temperature conditions and identified enzyme abundances that supported a high A TP yield. 159 \n  160 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\n 161 \n  162 \nFigure 3.  Changes in steady states of the modified model with A TP/ADP ratio homeostasis. (A) \nSteady-state concentration of D-glycerate 1,3-bisphosphate at each temperature. (B, C) Steady- state \nATP/ADP ratio (B) and A TP yield (C) at each temperature. (D) Relative changes in reaction rates at \nsteady state for each temperature. \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nA naïve approach of randomly varying enzyme levels and integrating the system dynamics 163 \nuntil convergence for each set would require computationally intensive simulations. To overcome this, 164 \nwe developed a faster method to sample numerous low-temperature steady states. Because every 165 \nreaction in our model is decomposed into elementary steps, each reaction rate function depends on at 166 \nmost one or two variables and follows mass-action kinetics. Consequently, the rate functions are at 167 \nmost quadratic polynomials, with the second variable, if present, corresponding to the concentration 168 \nof a small-molecule substrate. By fixing the concentrations of small molecules near their values in the 169 \noriginal steady state, the rate functions become linear, transforming the system of differential equations 170 \ninto a linear one. Steady states of this linear system can then be obtained simply by solving a linear 171 \nequation (equationௗ(7) in Materials and Methods). 172 \nThe linearized system has infinitely many steady states. To sample one, we first selected a 173 \nrandom reference point of enzyme abundances 𝑥௥௘௙  , near those of the original steady state. Each 174 \nconcentration in the original steady state was multiplied by a value drawn from a log-normal 175 \ndistribution over the interval (1/3, 3.0). Broadening this range reduces the fraction of sampled states 176 \nthat are stable, so we chose the widest interval that still yielded stable solutions with reasonable 177 \nfrequency. We then solved a linear optimization problem to identify the steady state closest to 𝑥௥௘௙ , 178 \nsubject to additional constraints: (1) the glucose transport reaction must proceed inward while other 179 \ntransport reactions act outward, and (2) the total enzyme amount must equal that of the original steady 180 \nstate. These transport constraints were imposed to ensure the biological plausibility of the sampled 181 \nsteady states. A mathematical formulation and detailed description of the sampling procedure are 182 \nprovided in Materials and Methods. This approach samples steady states more efficiently than the 183 \nnaïve method. Because stability was not enforced during sampling, we assessed it post hoc using the 184 \nJacobian matrix at each steady state. Materials and Methods also describe how steady states were 185 \nclassified as stable based on the Jacobian. 186 \nUsing this method, we sampled 5,000 steady states with rate constants corresponding to 0 ௗ°C and 187 \nvisualized the 94-dimensional enzyme profiles in two dimensions via principal component analysis 188 \n(PCA; Figureௗ4A). Among all computed solutions, the stable steady states were concentrated in a 189 \nspecific region, where TCA cycle enzyme levels were elevated. Figure ௗ4B shows the sum of TCA 190 \nenzyme abundances at each steady state in the PCA plot. These results suggest that increasing TCA 191 \ncycle enzyme concentrations is important for stabilizing the metabolic state under low-temperature 192 \nconditions. 193 \n  194 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\n 195 \n  196 \nFigure 4.  (A) PCA plot of 94- dimensional enzyme abundances at steady states obtained with \nparameters for 0 ௗ°C. Explained variance: PC1, 16.9%; PC2, 4.7%. Stability of each steady state is \nindicated by color: blue, stable; grey, unstable. (B) TCA cycle enzyme abundances (color- coded) at \nsteady states under low-temperature conditions. (C) Box plot of ATP yield for stable steady states at \neach temperature. \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nThe above results demonstrate that adjusting enzyme abundances allows the system to 197 \nachieve stable steady states near the original steady state, even under low-temperature conditions. 198 \nHowever, ATP yields at these stable states are not necessarily high. To evaluate A TP yield under low 199 \ntemperatures, we sampled stable steady states across temperatures from 37ௗ°C to 0ௗ°C using the same 200 \nmethod and assessed their A TP production. The A TP yield remained high even at low temperatures. 201 \nFigureௗ4C shows box plots of ATP yields for stable steady states at 37 ௗ°C, 19ௗ°C, and 0 ௗ°C. These 202 \nfindings indicate that cells can maintain a high ATP yield during a low-temperature shift by modulating 203 \nenzyme abundances. 204 \nIn practice, cells may not experience a sudden drop in temperature but rather undergo 205 \ngradual cooling. We therefore investigated which enzyme reallocations would allow the system to 206 \nmaintain a high ATP yield as temperature decreases from 37 °C to 0ௗ°C. To do this, we fixed the target 207 \nconcentrations of free enzymes and substrate-enzyme complexes ( 𝑥௥௘௙ ), and computed the steady-208 \nstate enzyme abundances at each temperature. Along this path, enzyme levels at stable steady states 209 \nchanged with temperature: TCA cycle enzymes increased, whereas ETC enzyme levels decreased. 210 \nFiguresௗ5A and 5B show the summed abundances of TCA cycle and ETC enzymes at stable steady 211 \nstates across temperatures. This trend is consistent with experimental observations that cold-adapted 212 \ncells upregulate TCA cycle enzymes and downregulate ATP synthase subunits (Knapp et ௗal., 2024; 213 \nPhadtare & Inouye, 2004). 214 \nBecause rate constants also vary with temperature, changes in enzyme abundance do not 215 \ntranslate directly into proportional changes in flux. Therefore, to understand how enzyme adjustments 216 \nmaintain a high ATP yield during cooling, it is essential to examine the resulting metabolic fluxes 217 \nrather than enzyme levels alone. We analyzed reaction rates normalized by glucose influx to 218 \ninvestigate flux-level modulations arising from enzyme redistribution. The analysis revealed that, as 219 \ntemperature decreased, normalized flux through the TCA cycle increased, along with fluxes of 220 \nreactions involved in ATP synthesis or consumption. Figuresௗ5C and 5D show the summed reaction 221 \nfluxes of the TCA cycle and ATP-utilizing reactions at stable steady states for each temperature. These 222 \nfindings indicate that the adjusted enzyme abundances enhanced ATP–ADP turnover, analogous to the 223 \nrapid ATP–ADP exchange artificially imposed in the previous section. Together, these results suggest 224 \nthat maintaining a high A TP/ADP ratio through increased cofactor turnover allows cells to sustain high 225 \nATP yield even at low temperatures. 226 \n  227 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\n  228 \nFigure 5. (A, B) Changes in the summed abundances of TCA cycle enzymes (A) and ETC enzymes \n(B) at stable steady states across temperatures. Each line represents one 𝒙𝒓𝒆𝒇. (C, D) Changes in the \nsummed absolute values of TCA cycle reaction rates (C) and A TP-utilizing reaction rates (D), both \nnormalized to glucose influx, at stable steady states across temperatures. Each line represents one \n𝒙𝒓𝒆𝒇. \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nDiscussion 229 \nIn this study, we computationally investigated how perturbations in environmental parameters 230 \ninfluence the stability of a metabolic state. By varying rate constants in a temperature-dependent 231 \nmanner, we found that lowering the temperature caused a bifurcation and a steep drop in the ATP yield. 232 \nThis finding provides a model-based illustration of how complex metabolic systems can be 233 \ndestabilized by parameter perturbations (Chakrabarti et al., 2013; Murayama et al., 2017; Steuer et al., 234 \n2006). 235 \nOne way to maintain a high ATP yield was to keep the ATP/ADP ratio constant through 236 \nexternal control. Several studies have reported that coenzymes such as ATP and ADP exert a major 237 \ninfluence on the perturbation responses of kinetic models (Hatakeyama & Furusawa, 2017; Himeoka 238 \n& Mitarai 2022 ; Himeoka & Furusawa, 2025). In our model, fixing this cofactor ratio preserved a 239 \nstable metabolic state across the temperature range and maintained a high A TP yield as shown in Figure 240 \n3C. We also demonstrated that the A TP yield can be maintained at a high level by adjusting enzyme 241 \nabundances. Under such enzyme-tuning scenarios, the relative flux through A TP-utilizing reactions 242 \nincreased, keeping the A TP/ADP ratio nearly constant. Thus, regulation of enzyme abundances can, in 243 \nprinciple, achieve the same effect as cofactor-ratio control. 244 \nAnalysis of enzyme abundances to maintain the ATP yield revealed an increase in the TCA 245 \ncycle enzymes and a decrease in ETC enzymes. This trend aligns with experimental reports that E. 246 \ncoli increases the TCA cycle enzyme levels and decreases ATP synthase subunits under cold 247 \ntemperature (Knapp et al., 2024; Phadtare & Inouye, 2004) . This suggests that cells may modulate 248 \nenzyme abundances in vivo  to maintain the ATP/ADP ratio and prevent an abrupt shift in their 249 \nmetabolic state. Nevertheless, these findings require experimental validation, such as comprehensive 250 \nfluxomic analyses that cover pathways beyond our model. 251 \nOur analysis perturbed only temperature, starting from parameters fitted to the 37ௗ°C metabolic state. 252 \nTherefore, the conclusions depend on that particular parameter set. Additionally, the model is limited 253 \nto central metabolism and omits lipid pathways and many amino-acid biosynthesis routes that are 254 \nreported to be important for cold adaptation (Phadtare & Inouye, 2008; Zhang & Gross, 2021) . 255 \nNetwork structure modifications in kinetic models are known to exert strong influences on responses 256 \nfor parameter perturbations (Hishida et al., 2023). Future work can probe cellular responses to cooling 257 \nwith higher quantitative accuracy using a more detailed network and more precise temperature 258 \ndependence of parameters.  259 \nThe same computational framework can be applied to environmental perturbations beyond 260 \ntemperature. By allowing kinetic parameters to vary with factors such as pH, osmotic pressure, or 261 \nnutrient availability, we can reveal which metabolic functions become unstable and what operations—262 \nsuch as cofactor balance control or targeted enzyme tuning—can restore stability. 263 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nOverall, our results highlight the latent risk of destabilization in complex metabolic systems 264 \nunder environmental perturbations and identify cofactor balance—specifically the ATP/ADP ratio—265 \nas a critical variable for maintaining stability. Applying this approach to other external parameters may 266 \nallow systematic identification of key variables that cells must regulate to preserve metabolic stability. 267 \nSuch insights would advance our understanding of the principles underlying biological robustness. 268 \n  269 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\n  270 \nSupplementary Figure S1.  Bifurcation diagram of small chemicals. The blue line represents the \nsteady state used for parameter fitting in Khodayari etௗal. (2014) at 37ௗ°C, while the red and grey lines \ncorrespond to steady states that exist below 14.2 ௗ°C. The steady state at the blue line shifts to that \nrepresented by the red line due to the temperature-induced bifurcation. \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\n 271 \nReaction Pathway Formula \n2DDA7Ptpp_ex   1 2dda7p <--->  \n3PGtpp_ex   1 3pg <--->  \nACALD \n \n 1 acald + 1 coa + 1 nad <---> 1 accoa + 1 nadh  \nACALDtpp_ex \n \n 1 acald =>  \nACK \n \n 1 ac + 1 atp <---> 1 actp + 1 adp  \nACONTa TCA cycle  1 cit <---> 1 acon-C + 1 h2o  \nACONTb TCA cycle  1 acon-C + 1 h2o <---> 1 icit  \nACS \n \n 1 ac + 1 atp + 1 coa <---> 1 accoa + 1 amp \nACt2rpp_ex_H \n \n 1 ac =>  \nADK \n \n 1 amp + 1 atp <---> 2 adp  \nAKGDH TCA cycle  1 akg + 1 coa + 1 nad <---> 1 nadh + 1 succoa  \nAKGt2rpp_ex \n \n 1 akg <--->  \nALCD \n \n 1 etoh + 1 nad <---> 1 acald + 1 nadh  \nATPM \n \n 1 atp <---> 1 adp \nATPS Electron \ntransport chain \n 1 adp <---> 1 atp \nCS TCA cycle  1 accoa + 1 oaa <---> 1 cit + 1 coa \nCYTBD Electron \ntransport chain \n 1 q8h2 <---> 1 q8  \nDDPA Pentose \nphosphate \npathway \n 1 e4p + 1 pep <---> 1 2dda7p \nDHAPtpp_ex \n \n 1 dhap  <---> \nD-LACt2pp_ex_H \n \n 1 lac-D =>  \nE4Ptpp_ex \n \n 1 e4p <---> \nEDA \n \n 1 2ddg6p <---> 1 g3p + 1 pyr \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nEDD \n \n 1 6pgc <---> 1 2ddg6p \nENO Glycolysis  1 2pg <---> 1 pep  \nETOHt2rpp_ex_H \n \n 1 etoh =>  \nF6Ptpp_ex \n \n 1 f6p <---> \nFBA Glycolysis  1 fdp <---> 1 dhap + 1 g3p  \nFBP Glycolysis  1 fdp <---> 1 f6p \nFBP_Glpx \n \n 1 fdp <---> 1 f6p  \nFORt2pp_ex_H \n \n 1 for =>  \nFORtppi_ex \n \n 1 for => \nFRD TCA cycle  1 fum + 1 q8h2 <---> 1 q8 + 1 succ  \nFRUpts2pp_ex \n \n 1 1 f6p + 1 pyr =>  1 pep  \nFUM TCA cycle  1 fum <---> 1 mal  \nFUMt2_2pp_ex_H \n \n 1 fum =>  \nG3Ptpp_ex \n \n 1 g3p <---> \nG6PDH Pentose \nphosphate \npathway \n 1 g6p + 1 nadp <---> 1 6pgl + 1 nadph  \nGAPD Glycolysis  1 g3p + 1 nad <---> 1 13dpg + 1 nadh  \nGLCptspp_ex_exi Glycolysis  1 pep <---> 1 g6p + 1 pyr  \nGLCS \n \n 1 adpglc <---> 1 adp+ 1 glycogen \nGLGC \n \n 1 atp + 1 g1p <---> 1 adpglc \nGLNS \n \n 1 atp + 1 glu <---> 1 adp + 1 gln  \nGLUD \n \n 1 glu + 1 nadp <---> 1 akg + 1 nadph  \nGLUN \n \n 1 gln <---> 1 glu  \nGLUS \n \n 1 akg + 1 gln + 1 nadph <---> 2 glu + 1 nadp  \nGLUt2rpp_ex_H \n \n 1 glu =>  \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nGND Pentose \nphosphate \npathway \n 1 6pgc + 1 nadp <---> 1 nadph + 1 ru5p  \nHEX \n \n 1 atp + 1 glc-D <---> 1 adp + 1 g6p  \nICDH TCA cycle  1 icit + 1 nadp <---> 1 akg + 1 nadph  \nICL \n \n 1 icit <---> 1 glx + 1 succ  \nLDH \n \n 1 lac-D + 1 nad <---> 1 nadh + 1 pyr  \nMALS \n \n 1 accoa + 1 glx <---> 1 coa + 1 mal  \nMALt2_2pp_ex_H \n \n 1 mal =>  \nMDH TCA cycle  1 mal + 1 nad <---> 1 nadh + 1 oaa  \nMDH2 TCA cycle  1 mal + q8 <---> oaa + q8h2 \nME1 \n \n 1 mal + 1 nad <---> 1 nadh + 1 pyr  \nME2 \n \n 1 mal + 1 nadp <---> 1 nadph + 1 pyr  \nNADH Electron \ntransport chain \n 1 nadh + 1 q8 <---> 1 nad + 1 q8h2  \nNADTRHD \n \n 1 nad + 1 nadph <---> 1 nadh + 1 nadp  \nOAAtpp_ex \n \n 1 oaa <---> \nPDH Glycolysis  1 coa + 1 nad + 1 pyr <---> 1 accoa + 1 nadh  \nPEPtpp_ex \n \n 1 pep <---> \nPFKA Glycolysis  1 atp + 1 f6p <---> 1 adp + 1 fdp  \nPFKB Glycolysis  1 atp + 1 f6p <---> 1 adp + 1 fdp \nPFL \n \n 1 coa + 1 pyr <---> 1 accoa + 1 for  \nPGI Glycolysis  1 g6p <---> 1 f6p  \nPGK Glycolysis  1 13dpg + 1 adp <---> 1 3pg + 1 atp \nPGL Pentose \nphosphate \npathway \n 1 6pgl <---> 1 6pgc  \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nPGL_spon Pentose \nphosphate \npathway \n 1 6pgl <---> 1 6pgc  \nPGM Glycolysis  1 2pg <---> 1 3pg  \nPGMT \n \n 1 g1p <--->  1 g6p \nPPC \n \n 1 pep <---> 1 oaa  \nPPCK \n \n 1 atp + 1 oaa <---> 1 adp + 1 pep  \nPPS Glycolysis  1 atp + 1 pyr <---> 1 amp + 1 pep  \nPTA \n \n 1 accoa <---> 1 actp + 1 coa  \nPYKA Glycolysis  1 adp + 1 pep <---> 1 atp + 1 pyr \nPYKF Glycolysis  1 adp + 1 pep <---> 1 atp + 1 pyr  \nPYRpp_ex \n \n 1 pyr <---> \nPYRt2rpp_ex_H \n \n 1 pyr <---> \nR5PP Pentose \nphosphate \npathway \n 1 r5p <---> 1 rib \nRIB_Dtpp_ex \n \n 1 rib <---> \nRPE Pentose \nphosphate \npathway \n 1 ru5p <---> 1 xu5p  \nRPI Pentose \nphosphate \npathway \n 1 r5p <---> 1 ru5p  \nRU5Ptpp_ex \n \n 1 ru5p => \nS7Ptpp_ex \n \n 1 s7p => \nSUCCt2_2pp_ex_H \n \n 1 succ => \nSUCCt3pp_ex_H \n \n 1 succ <---> \nSUCD TCA cycle  1 q8 + 1 succ <---> 1 fum + 1 q8h2  \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nSUCOAS TCA cycle  1 atp + 1 coa + 1 succ <---> 1 adp + 1 succoa  \nTALA Pentose \nphosphate \npathway \n 1 g3p + 1 s7p <---> 1 e4p + 1 f6p  \nTHD \n \n 1 nadh + 1 nadp <---> 1 nad + 1 nadph  \nTKT1 Pentose \nphosphate \npathway \n 1 r5p + 1 xu5p <---> 1 g3p + 1 s7p  \nTKT2 Pentose \nphosphate \npathway \n 1 e4p + 1 xu5p <---> 1 f6p + 1 g3p  \nTPI Glycolysis  1 dhap <---> 1 g3p  \n 272 \nSupplementary Table S1. List of all reactions in the model used in our simulations. 273 \n  274 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nReaction Regulator Regulation Type \nACALD amp Inhibition \nEDA 6pgc Inhibition \nEDA g3p Inhibition \nFBA dhap Inhibition \nFBP adp Inhibition \nFBP amp Inhibition \nFBP pi Inhibition \nFBP pep Inhibition \nFUM cit Inhibition \nG6PDH nadh Inhibition \nG6PDH nadph Inhibition \nGLGC amp Inhibition \nGLUN adp Inhibition \nGLUN atp Inhibition \nGLUN glu Inhibition \nGLUN nad Inhibition \nGLUN nadh Inhibition \nGLUN nadp Inhibition \nGLUN nadph Inhibition \nGLUN pi Inhibition \nGLUS glu Inhibition \nGND atp Inhibition \nGND ru5p Inhibition \nGND fdp Inhibition \nGND nadp Inhibition \nICDH glx Inhibition \nICDH oxa Inhibition \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nICDH pep Inhibition \nICL akg Inhibition \nICL acon Inhibition \nICL pep Inhibition \nICL succ Inhibition \nLDH atp Inhibition \nPGI 6pgc Inhibition \nPGI pep Inhibition \nPGMT accoa Inhibition \nPGMT g1p Inhibition \nPPC mal Inhibition \nPPC pep Inhibition \nPPCK atp Inhibition \nPPCK pep Inhibition \nPPS akg Inhibition \nPPS adp Inhibition \nPPS amp Inhibition \nPPS oaa Inhibition \nPPS pep Inhibition \nPTA adp Inhibition \nPTA atp Inhibition \nPTA nadh Inhibition \nPTA nadph Inhibition \nPTA pep Inhibition \nPTA pyr Inhibition \nPYKA PYKF atp Inhibition \nPYKA PYKF succoa Inhibition \nRPI amp Inhibition \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nTALA pi Inhibition \n 275 \nSupplementary Table S2. List of all regulations in the model used in our simulations. 276 \n 277 \n 278 \n  279 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nMaterials and Methods 280 \n 281 \nModel Description 282 \nThe model used in this analysis is based on the model from Himeoka et al. (2025), which is a modified 283 \nversion of the Khodayari et al. (2014) model, adjusted to ensure the existence of a stable steady state. 284 \nIn the present study, the parameters of inflow and outflow reactions were further modified so that 285 \nglucose serves as the sole carbon source in the temperature range used for the simulations. The model 286 \nincludes glycolysis, the TCA cycle, the pentose phosphate pathway (PPP), glutamine synthesis, the 287 \nelectron transport chain (ETC), and the Entner-Doudoroff pathway. See Supplementary Table S1 and 288 \nS2 for complete reaction formulae, pathway components, and regulations. 289 \nIn our model, all reactions are decomposed into elementary reactions that represent 290 \nsubstrate-enzyme binding and dissociation. This is illustrated using a simple reversible reaction as 291 \nfollows: 292 \nFor a reversible reaction A + B ⇌ C catalyzed by an enzyme E, we write: 293 \nA + E ⇌ EA 294 \nEA + B ⇌ EAB 295 \nEAB ⇌ EC 296 \nEC ⇌ E + C. 297 \nAll steps are reversible, and the forward and reverse directions are treated as distinct elementary 298 \nreactions. The total concentration of free and complex enzymes does not change in time: 299 \n[𝐸] + [𝐸𝐴] + [𝐸𝐴𝐵] + [𝐸𝐶] = const. 303 \nThe rates of these elementary reactions are calculated based on mass-action kinetics. Note that in our 300 \nmodel, chemical concentrations are dimensionless because they are normalized to 1.0 at the reference 301 \nsteady state as a standard procedure in ensemble modeling (Khodayari et al., 2014; Tran et al., 2008). 302 \nIn our model, the time evolution of concentrations is described by an ODE system 304 \n𝑑𝑥(𝑡)\n𝑑𝑡 = 𝑆𝑟൫𝑥(𝑡); 𝑘(𝑇)൯, (1) 311 \nwhere 𝑥  is the concentration vector of small molecules, free enzymes, and substrate-enzyme 305 \ncomplexes; 𝑆 is the stoichiometric matrix; 𝑘(𝑇) is the set of rate constants at temperature 𝑇; and 306 \n𝑟(𝑥) is the vector of reaction rate functions. All reactions are decomposed into elementary reactions, 307 \nand each reaction rate is modeled by the law of mass action  308 \n𝑟௜(𝑥) = 𝑘௜(𝑇) ෑ 𝑥௝\nఔ ೕ೔\n௡\n௝ୀଵ\n, (2) 312 \nwith exponents 𝜈௝௜ ∈ {0,1} determined by the stoichiometry and 𝑘௜(𝑇) the temperature-dependent 309 \nrate constants (specified in the next section).  310 \nODE simulations were performed using the ode15s function in MATLAB (The MathWorks 313 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nInc., 2023). 314 \n 315 \nTemperature-Dependent Rate Constants 316 \nRate constants were adjusted with the Eyring-Polanyi equation 317 \n𝑘(𝑇) = 𝑘஻ 𝑇\nℎ exp ቆ − ∆𝐺‡\n𝑅𝑇 ቇ  (3) 324 \nfor each temperature. Here, 𝑘 is the rate constant, ℎ is the Planck constant, 𝑘஻  is the Boltzmann 318 \nconstant, 𝑅 is the gas constant, and 𝑇 is the absolute temperature (𝐾). The activation free energy 319 \n∆𝐺‡ for each reaction is back-calculated from parameter values obtained in Khodayari et al. 2014 320 \nusing the following equation: 321 \n∆𝐺‡ = −𝑅𝑇௥௘௙ ln ቆ 𝑘∗ ℎ\n𝑘஻ 𝑇௥௘௙\nቇ , (4) 325 \nwith 𝑇௥௘௙ = 310.15, and 𝑘∗ is the rate constant at 𝑇௥௘௙ . Based on this value, reaction rate constants 322 \nat altered temperatures were calculated. 323 \n 326 \nMaintaining the ATP/ADP Ratio 327 \nTo forcibly maintain the ATP/ADP ratio while lowering the temperature, we added a new reversible 328 \nreaction that exchanges A TP and ADP to the model. The rate constant for this exchange was set to 329 \n10଺ , a value larger than the median of all rate constants at 37 °C. During low-temperature simulations, 330 \nthis constant was not adjusted for temperature dependence. The concentrations of the free metabolites 331 \nat the steady state are normalized to unity in the ensemble modeling procedure carried out for the 332 \nmodel construction in Khodayari et al. (2014). Thus, ATP/ADP ratio in the original model is 1:1 in the 333 \nscale of normalized concentrations. Therefore, the force-balancing A TP/ADP ratio to unity 334 \ncorresponds to the recovery of the ATP/ADP ratio in the original model.  335 \n 336 \nPerturbed Initial Condition Sampling Procedure 337 \nTo identify multiple stable steady states in the system, we perform numerical integrations starting from 338 \nperturbed initial conditions with specified temperature parameters. These initial conditions are 339 \ngenerated by multiplying the concentrations of an original steady state by random values with ensuring 340 \nthat the conserved quantities are not changed. 341 \nIf a vector 𝑑 satisfies 342 \n𝑑୘ 𝑆 = 0, 345 \nthen the inner product of 𝑑 and a concentration vector 𝑥 does not change in time, since  343 \n𝑑\n𝑑𝑡 (𝑑୘ 𝑥) = 0. 346 \n 344 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nIf random perturbations are applied directly to the concentrations of the original steady state 347 \n𝑥௢௥௚ , the resulting concentration vector 𝑥௥௔௡ௗ  will generally have conserved quantities different from 348 \nthose of 𝑥௢௥௚ . To obtain perturbed concentrations that preserve the conserved quantities of 𝑥௢௥௚ , we 349 \nsearch for the point closest to 𝑥௥௔௡ௗ  that satisfies the conservation constraints.  350 \nThis is formulated as the following optimization problem: 351 \n𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 ‖𝑥 − 𝑥௥௔௡ௗ ‖ଵ  354 \n𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 355 \n𝑑௜\n୘ 𝑥 = 𝑑௜\n୘ 𝑥௢௥௚ , for all 𝑖. 356 \nWhere {𝑑௜}  form a basis of ker 𝑆୘ .  This optimization is done using Gurobi Optimizer (Gurobi 352 \nOptimization, 2024). 353 \n 357 \nSteady State Sampling Procedure 358 \nTo identify enzyme abundances that yield a stable metabolic state with low-temperature parameters, 359 \nwe sampled steady states for each temperature-specific parameter set and computed the corresponding 360 \nenzyme abundances. Each reaction is decomposed into elementary reactions whose reactants are (i) 361 \nsubstrate and enzyme, (ii) substrate and substrate-enzyme complex, or (iii) substrate-enzyme complex 362 \nonly. Therefore, each reaction rate function is a first- or second-order monomial in its variables. 363 \nFixing concentrations of substrates (small molecules) at their values near the original steady 364 \nstate turns every rate into a linear function of the free-enzyme or substrate-enzyme complex 365 \nconcentrations. 𝑥௘௡௭௬௠௘   is a vector of the remaining variables, i.e., the concentrations of free 366 \nenzymes and substrate-enzyme complexes. Using a matrix 𝑉, whose non-zero values are determined 367 \nby the temperature-specific rate constants and the fixed substrate concentrations, we have  368 \n𝑟 = 𝑉𝑥௘௡௭௬௠௘ (5) 370 \nwhere 𝑟 is the vector of reaction rates.  369 \nFor example, if the system has only one reaction A + B ⇌ C catalyzed by an enzyme E, 371 \nreaction rate functions are 372 \n𝑟 =\n⎝\n⎜\n⎜\n⎜\n⎜\n⎜\n⎛\n𝑘ଵ௙ [𝐴][𝐸]\n𝑘ଵ௕ [𝐸𝐴]\n𝑘ଶ௙ [𝐸𝐴][𝐵]\n𝑘ଶ௕ [𝐸𝐴𝐵]\n𝑘ଷ௙ [𝐸𝐴𝐵]\n𝑘ଷ௕ [𝐸𝐶]\n𝑘ସ௙ [𝐸𝐶]\n𝑘ସ௕ [𝐸][𝐶] ⎠\n⎟\n⎟\n⎟\n⎟\n⎟\n⎞\n, 375 \nwhere 𝑘௜ is rate constant of forward or backward direction of each elementary reaction. When the 373 \nsubstrate concentrations are fixed at ([𝐴]∗, [𝐵]∗, [𝐶]∗)୘ ,  the rate functions can be rewritten as 374 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\n𝑟 =\n⎝\n⎜\n⎜\n⎜\n⎜\n⎜\n⎛\n𝑘ଵ௙ [𝐴]∗ 0 0 0\n0 𝑘ଵ௕ 0 0\n0 𝑘ଶ௙ [𝐵]∗ 0 0\n0 0 𝑘ଶ௕ 0\n0 0 𝑘ଷ௙ 0\n0 0 0 𝑘ଷ௕\n0 0 0 𝑘ସ௙\n𝑘ସ௕ [𝐶]∗ 0 0 0 ⎠\n⎟\n⎟\n⎟\n⎟\n⎟\n⎞\n൮\n[𝐸]\n[𝐸𝐴]\n[𝐸𝐴𝐵]\n[𝐸𝐶]\n൲ . 376 \nUsing submatrix of the stoichiometric matrix corresponding to the remaining variables, 377 \ndenoted as 𝑆௘௡௭௬௠௘ , the time-course change of 𝑥௘௡௭௬௠௘  can be written as  378 \n𝑑𝑥௘௡௭௬௠௘\n𝑑𝑡 = 𝑆௘௡௭௬௠௘ 𝑉𝑥௘௡௭௬௠௘ . (6) 380 \nTherefore, steady-state solutions are obtained just by solving  379 \n𝑆௘௡௭௬௠௘ 𝑉𝑥௘௡௭௬௠௘ = 0. (7) 381 \nEq. (7) is less-determined linear equation when 𝑆௘௡௭௬௠௘ 𝑉 is not full-rank. In that case, Eq. 382 \n(7) has infinitely many solutions. Among such solutions of Eq. (7), we select the one closest to a 383 \nreference point ൫𝑥௥௘௙ ൯, which is randomly sampled around the enzyme abundances of the original 384 \nsteady state. We also require steady states that  385 \n1. Glucose is a single energy resource, and transport reactions of other chemicals work as efflux. 386 \n2. The total enzyme abundance of all reactions is equal to that of the original steady state.  387 \nSteady states satisfying those conditions are obtained by solving the following optimization problem:  388 \n𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 ‖𝑥௘௡௭௬௠௘ − 𝑥௥௘௙ ‖ଵ  392 \n𝑆𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 393 \n𝑆௘௡௭௬௠௘ 𝑉𝑥௘௡௭௬௠௘ = 0 394 \n෍ 𝑥௘௡௭௬௠௘ ,௜\n௜\n= const. 395 \n𝑣ீ௟௨௖௢௦௘ூ௡௣௨௧ > 0 396 \n𝑣ை௧௛௘௥்௥௔௡௦௣ < 0 397 \nThe solution 𝑥௘௡௭௬௠௘   obtained from this optimization, together with the fixed small molecule 389 \nconcentrations near the 37 °C steady state, yields a steady state of the full system. These optimizations 390 \nare done using Gurobi Optimizer (Gurobi Optimization, 2024). 391 \n 398 \nStability Assessment of Steady States 399 \nThe stability of each steady state found by the optimization was determined by computing the 400 \neigenvalues of the system’s Jacobi matrix at the steady state. The model contains 98 conserved 401 \nquantities, so the full 643 × 643 Jacobi matrix has 98 zero eigenvalues. The zero eigenvalues are often 402 \ntricky to distinguish from the very small, non-zero eigenvalues due to numerical errors. Thus, we 403 \navoided obtaining zero eigenvalues by pre-solving the dependent variables using the conserved 404 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nquantities. Note that one conserved quantity leads to a single linear equation. For each conserved 405 \nquantity, we choose a single chemical species and solve for its concentration as a function of the 406 \nremaining species’ concentrations. From this pre-solved ODE system, we obtain the reduced 545 × 407 \n545 Jacobi matrix. A steady state was classified as unstable if the largest real part of any eigenvalue 408 \nof the reduced Jacobi matrix was positive. All eigenvalues were computed with interval-arithmetic 409 \nroutines in INTLAB, ensuring that numerical errors did not flip the sign of any eigenvalue (Rump, 410 \n1999). 411 \n 412 \nAcknowledgements 413 \nThis work is supported in part by JSPS KAKENHI (23KJ1324 to A . H.; 22K15069, 24H01118 and 414 \n25H01390 to Y . H.; 22K21344 and 23H0247 to C. F.), and GteX Program Japan Grant No. 415 \nJPMJGX23B4. 416 \n  417 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nReferences 418 \nArcus, V . L., & Mulholland, A. J. (2020). Temperature, Dynamics, and Enzyme-Catalyzed Reaction 419 \nRates. Annual Review of Biophysics, 49(1), 163–180. https://doi.org/10.1146/annurev-biophys-420 \n121219-081520 421 \nChakrabarti, A., Miskovic, L., Soh, K. C., & Hatzimanikatis, V . (2013). Towards kinetic modeling of 422 \ngenome-scale metabolic networks without sacrificing stoichiometric, thermodynamic and 423 \nphysiological constraints. Biotechnology Journal , 8(9), 1043–1057. 424 \nhttps://doi.org/10.1002/biot.201300091 425 \nChristodoulou, D., Link, H., Fuhrer, T., Kochanowski, K., Gerosa, L., & Sauer, U. (2018). Reserve 426 \nFlux Capacity in the Pentose Phosphate Pathway Enables Escherichia coli’s Rapid Response to 427 \nOxidative Stress. Cell Systems, 6(5), 569-578.e7. https://doi.org/10.1016/j.cels.2018.04.009 428 \nEl-Samad, H., Kurata, H., Doyle, J. C., Gross, C. A., & Khammash, M. (2005). Surviving heat shock: 429 \nControl strategies for robustness and performance. Proceedings of the National Academy of 430 \nSciences, 102(8), 2736–2741. https://doi.org/10.1073/pnas.0403510102 431 \nEvans, M. G., & Polanyi, M. (1935). Some applications of the transition state method to the calculation 432 \nof reaction velocities, especially in solution. Transactions of the Faraday Society , 31, 875. 433 \nhttps://doi.org/10.1039/tf9353100875 434 \nEyring, H. (1935). The Activated Complex in Chemical Reactions. The Journal of Chemical Physics, 435 \n3(2), 107–115. https://doi.org/10.1063/1.1749604 436 \nGruber, C. H., Noor, E., Buffing, M. F., & Sauer, U. (2025). Systematic identification of allosteric 437 \neffectors in Escherichia coli metabolism. Proceedings of the National Academy of Sciences of 438 \nthe United States of America, 122(10). https://doi.org/10.1073/pnas.2423767122 439 \nGurobi Optimization, L. (2024). Gurobi Optimizer Reference Manual. 440 \nHatakeyama, T. S., & Furusawa, C. (2017). Metabolic dynamics restricted by conserved carriers: 441 \nJamming and feedback. PLoS Computational Biology , 13(11). 442 \nhttps://doi.org/10.1371/journal.pcbi.1005847 443 \nHimeoka, Y ., & Mitarai, N. (2022). Emergence of Growth and Dormancy from a Kinetic Model of the 444 \nEscherichia Coli Central Carbon Metabolism., Physical Review Research 4 (4): 043223. 445 \nHimeoka, Y ., & Furusawa, C. (2025). Perturbation-response analysis of in silico metabolic dynamics 446 \nrevealed hard-coded responsiveness in the cofactors and network sparsity. ELife, 13. 447 \nhttps://doi.org/10.7554/eLife.98800 448 \nHishida, A., Okada, T., & Mochizuki, A. (2023). Patterns of change in regulatory modules of chemical 449 \nreaction systems induced by network modification. PNAS Nexus , 3(1). 450 \nhttps://doi.org/10.1093/pnasnexus/pgad441 451 \nJozefczuk, S., Klie, S., Catchpole, G., Szymanski, J., Cuadros-Inostroza, A., Steinhauser, D., Selbig, 452 \nJ., & Willmitzer, L. (2010). Metabolomic and transcriptomic stress response of Escherichia coli. 453 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nMolecular Systems Biology, 6. https://doi.org/10.1038/msb.2010.18 454 \nKhodayari, A., Zomorrodi, A. R., Liao, J. C., & Maranas, C. D. (2014). A kinetic model of Escherichia 455 \ncoli core metabolism satisfying multiple sets of mutant flux data. Metabolic Engineering, 25, 456 \n50–62. https://doi.org/10.1016/j.ymben.2014.05.014 457 \nKnapp, B. D., Willis, L., Gonzalez, C., Vashistha, H., Jammal-Touma, J., Tikhonov, M., Ram, J., 458 \nSalman, H., Elias, J. E., & Huang, K. C. (2024). Metabolic rearrangement enables adaptation of 459 \nmicrobial growth rate to temperature shifts. Nature Microbiology , 10(1), 185–201. 460 \nhttps://doi.org/10.1038/s41564-024-01841-4 461 \nMiyakoshi, M. (2024). Multilayered regulation of amino acid metabolism in Escherichia coli. Current 462 \nOpinion in Microbiology, 77, 102406. https://doi.org/10.1016/j.mib.2023.102406 463 \nMurayama, Y ., Kori, H., Oshima, C., Kondo, T., Iwasaki, H., & Ito, H. (2017). Low temperature 464 \nnullifies the circadian clock in cyanobacteria through Hopf bifurcation. Proceedings of the 465 \nNational Academy of Sciences of the United States of America , 114(22), 5641–5646. 466 \nhttps://doi.org/10.1073/pnas.1620378114 467 \nNoll, P., Lilge, L., Hausmann, R., & Henkel, M. (2020). Modeling and Exploiting Microbial 468 \nTemperature Response. Processes, 8(1), 121. https://doi.org/10.3390/pr8010121 469 \nPhadtare, S. (2004). Recent Developments in Bacterial Cold-Shock Response. Current Issues in 470 \nMolecular Biology. https://doi.org/10.21775/cimb.006.125 471 \nPhadtare, S., & Inouye, M. (2004). Genome-wide transcriptional analysis of the cold shock response 472 \nin wild-type and cold-sensitive, quadruple-csp-deletion strains of Escherichia coli. Journal of 473 \nBacteriology, 186(20), 7007–7014. https://doi.org/10.1128/JB.186.20.7007-7014.2004 474 \nPhadtare, S., & Inouye, M. (2008). The Cold Shock Response. EcoSal Plus , 3(1). 475 \nhttps://doi.org/10.1128/ecosalplus.5.4.2 476 \nRiccardi, C., Calvanese, M., Ghini, V ., Alonso-Vásquez, T., Perrin, E., Turano, P., Giurato, G., Weisz, 477 \nA., Parrilli, E., Tutino, M. L., & Fondi, M. (2023). Metabolic Robustness to Growth Temperature 478 \nof a Cold- Adapted Marine Bacterium. MSystems, 8(2). 479 \nhttps://doi.org/10.1128/msystems.01124-22 480 \nRump, S. M. (1999). INTLAB — INTerval LABoratory. In Developments in Reliable Computing (pp. 481 \n77–104). Springer Netherlands. https://doi.org/10.1007/978-94-017-1247-7_7 482 \nSawers, R. G. (2025). How FocA facilitates fermentation and respiration of formate by Escherichia 483 \ncoli. In Journal of Bacteriology  (V ol. 207, Issue 2). American Society for Microbiology. 484 \nhttps://doi.org/10.1128/jb.00502-24 485 \nSchmidt, A., Kochanowski, K., Vedelaar, S., Ahrné, E., V olkmer, B., Callipo, L., Knoops, K., Bauer, 486 \nM., Aebersold, R., & Heinemann, M. (2016). The quantitative and condition-dependent 487 \nEscherichia coli proteome. Nature Biotechnology , 34(1), 104–110. 488 \nhttps://doi.org/10.1038/nbt.3418 489 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint \n\nStephanopoulos, G. ., Aristidou, A. A. ., & Nielsen, J. Høiriis. (1998). Metabolic engineeringࣟ :490 \nprinciples and methodologies. Academic Press. 491 \nSteuer, R., Gross, T., Selbig, J., & Blasius, B. (2006). Structural kinetic modeling of metabolic 492 \nnetworks. Proceedings of the National Academy of Sciences , 103(32), 11868–11873. 493 \nhttps://doi.org/10.1073/pnas.0600013103 494 \nStrutz, J., Martin, J., Greene, J., Broadbelt, L., & Tyo, K. (2019). Metabolic kinetic modeling provides 495 \ninsight into complex biological questions, but hurdles remain. Current Opinion in Biotechnology, 496 \n59, 24–30. https://doi.org/10.1016/j.copbio.2019.02.005 497 \nSudarsan, S., Demling, P., Ozdemir, E., Ben Ammar, A., Mennicken, P., Buescher, J. M., Meurer, G., 498 \nEbert, B. E., & Blank, L. M. (2025). Acetol biosynthesis enables NADPH balance during 499 \nnitrogen limitation in engineered Escherichia coli. Microbial Cell Factories , 24(1). 500 \nhttps://doi.org/10.1186/s12934-025-02687-z 501 \nThe MathWorks Inc. (2023). MATLAB version: R2023b. The MathWorks Inc. 502 \nTran, L. M., Rizk, M. L., & Liao, J. C. (2008). Ensemble modeling of metabolic networks. Biophysical 503 \nJournal, 95(12), 5606–5617. https://doi.org/10.1529/biophysj.108.135442 504 \nWendering, P., & Nikoloski, Z. (2023). Model-driven insights into the effects of temperature on 505 \nmetabolism. Biotechnology Advances , 67, 108203. 506 \nhttps://doi.org/10.1016/j.biotechadv.2023.108203 507 \nYe, Y ., Zhang, L., Hao, F., Zhang, J., Wang, Y ., & Tang, H. (2012). Global metabolomic responses of 508 \nescherichia coli to heat stress. Journal of Proteome Research , 11(4), 2559–2566. 509 \nhttps://doi.org/10.1021/pr3000128 510 \nZhang, Y ., & Gross, C. A. (2021). Cold Shock Response in Bacteria. Annual Review of Genetics, 55(1), 511 \n377–400. https://doi.org/10.1146/annurev-genet-071819-031654 512 \nZhao, J., Chen, K., Palsson, B. O., & Yang, L. (2024). StressME: Unified computing framework of 513 \nEscherichia coli metabolism, gene expression, and stress responses. PLoS Computational 514 \nBiology, 20(2 February). https://doi.org/10.1371/journal.pcbi.1011865 515 \n  516 \n 517 \n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}