Analysis of Metabolic Stability under Environmental Perturbations Using a Kinetic Model

13 Cells robustly maintain metabolic function despite environmental fluctuations affecting reaction 14 kinetics, yet kinetic models often exhibit fragility. We investigated this discrepancy by examining the 15 effects of temperature changes on a kinetic model of Escherichia coli central metabolism. A gradual 16 temperature decrease caused an abrupt transition to a state with sharply reduced ATP production 17 efficiency, triggered by an elevated ATP/ADP ratio creating a glycolysis bottleneck. Introducing a 18 rapid ATP-ADP exchange reaction prevented this ATP/ADP ratio surge and sustained high ATP 19 production efficiency across temperatures. Furthermore, altering enzyme abundances similarly 20 maintained high A TP production, with predicted enzyme regulation aligning with experimental 21 observations in E. coli at low temperatures. Our findings indicate that balancing key cofactors, such 22 as the ATP/ADP ratio, is crucial for preserving metabolic stability under environmental perturbations. 23 24 25 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint

26 Cellular metabolism comprises thousands of reactions that together generate energy, build 27 biomolecules, and recycle cellular components. As external conditions fluctuate, cells adjust their 28 metabolic state to maintain function. For example, bacteria switch between fermentative and 29 respiratory regimes, reroute carbon through the pentose phosphate pathway (PPP) under oxidative 30 stress, and upregulate anaplerotic reactions when nutrients are scarce (Christodoulou et al., 2018; 31 Sawers, 2025; Sudarsan et al., 2025). This plasticity arises from regulatory mechanisms such as 32 enzyme reallocation and allosteric modulation, enabling cells to occupy distinct, environment-specific 33 metabolic states while sustaining growth (Gruber et al., 2025; Miyakoshi, 2024). 34 Mathematical models of metabolic dynamics have been widely used to probe this plasticity 35 (El-Samad et al., 2005; Stephanopoulos et al., 1998; Strutz et al., 2019). Such models allow analysis 36 of how steady states shift as parameters vary. A recurring finding is that steady-state stability is often 37 fragile to parameter changes. For instance, a kinetic model of E. coli core metabolism was stable for 38 only a small fraction of flux and concentration sets consistent with its kinetic parameters (Chakrabarti 39 et al., 2013). Similarly, Steuer et al. (2006) sampled thousands of parameter sets for a kinetic model 40 of yeast glycolysis and found that over 80% of resulting steady states were unstable. These narrow 41 stability ranges suggest that environmental perturbations can readily destabilize cellular metabolism. 42 Yet, despite the fragility predicted by models, real cells maintain stable metabolic states as 43 temperature and nutrient availability change (Phadtare, 2004; Phadtare & Inouye, 2008; Riccardi et 44 al., 2023; Zhang & Gross, 2021). This discrepancy implies the existence of mechanisms that stabilize 45 metabolism. Indeed, cells adjust enzyme levels and pathway activities in response to environmental 46 cues (Jozefczuk et al., 2010; Schmidt et al., 2016; Ye et al., 2012; Zhao et al., 2024). However, how 47 such regulation preserves dynamic stability remains unclear. 48 Among environmental perturbations, temperature shifts directly and broadly affect 49 metabolism (Noll et al., 2020; Wendering & Nikoloski, 2023). Temperature modulates the rate 50 constants of all intracellular reactions, with magnitudes set by each reaction’s activation free energy 51 (Evans & Polanyi, 1935; Eyring, 1935). Because forward and reverse reactions need not respond 52 identically, temperature changes can alter both rates and reversibility, reshaping flux throughout the 53 network. Such nonuniform effects may destabilize steady states. Clarifying how cells prevent this 54 outcome would reveal general principles of metabolic robustness. 55 Here, we investigate mechanisms that keep complex metabolic systems stable under 56 temperature-induced parameter perturbations. Using a kinetic model of E. coli central metabolism 57 constructed from experimental data, we varied model parameters continuously with temperature and 58 observed a bifurcation: the system transitioned discontinuously to a new steady state. The post-59 bifurcation state exhibited an elevated A TP/ADP ratio and reduced A TP production efficiency. 60 Imposing operations that maintain the ATP/ADP ratio during the temperature shift eliminated the 61 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint bifurcation and preserved ATP production efficiency even at low temperatures. Moreover, adjusting 62 enzyme abundances to increase A TP-utilizing reaction rates stabilized the ATP/ADP ratio. These 63

indicate that maintaining the balance of key cofactors is critical for stabilizing metabolic 64 systems under parameter perturbations caused by environmental change. 65 66

67 Temperature decrease induces a bifurcation 68 To examine how a low-temperature shift affects the metabolic state of E. coli central metabolism, we 69 simulated metabolic dynamics using rate constants adjusted for different temperatures. Figure 1 shows 70 the reaction network of the kinetic model used in this study (see Materials and Methods for details). 71 Rate constants at each temperature were calculated using the Eyring–Polanyi equation. We focused on 72 temperature downshifts because heating can cause irreversible enzyme inactivation, which is not 73 captured by the standard Eyring–Polanyi formulation (Arcus & Mulholland, 2020; Noll et al., 2020). 74 Activation free energies (∆𝐺‡) for individual reactions were obtained from parameter values fitted to 75 the in vivo metabolic state at 37 °C reported by Khodayari et al. (2014) (see Materials and Methods 76 for the details). Without a temperature shift, the system converges to a stable steady state in which 77 glucose serves as the sole carbon source and ATP is produced via glycolysis, the TCA cycle, and the 78 electron transport chain (ETC). We refer to this as the “original steady state” throughout the paper. 79 When the temperature is lowered from 37 °C to 0 °C, each rate constant decreases according 80 to its activation free energy. As a result, the distribution of rate constants shifts, and the median value 81 decreases by about 10-fold (Figure 2A). We simulated system dynamics using the rate constants at 82 0 °C, starting from the original steady state. Numerical simulations showed that metabolite 83 concentrations converged to a new stable steady state (Figure 2B). Compared with the original state, 84 metabolite concentrations in this new steady state changed by factors ranging from 10ିଽ to 10ସ —85 far greater than the fold change observed in the median rate constant. 86 To determine the temperature at which this substantial shift occurs, we tracked the steady 87 state as the temperature was gradually decreased from 37 °C to 0 °C. At each step, system dynamics 88 were integrated up to 𝑡 = 10ଵ଴ 𝑠 to ensure convergence, and the resulting concentrations were used 89 as initial conditions for the next step, with temperature reduced by 0.1 °C. At 14.2 °C, the steady state 90 shifted discontinuously, with abrupt changes in the concentrations of all metabolites (Figure 2C shows 91 one example; Supplementary Figure S1 presents all metabolite concentrations at the converged state 92 across temperatures). Reaction rates also changed discontinuously, and relative changes compared 93 with the original steady state are shown in Figure 2D. 94 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint 95 96 Figure 1. Reaction network of the kinetic model used in our simulations. Although each reaction was decomposed into elementary steps for the simulations, the network is shown here with the pre - decomposition reactions for clarity. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint This discontinuous shift could arise from either the disappearance or destabilization of the 97 original steady state, or from a switch between coexisting attractors. To discriminate between these 98 possibilities, we enumerated multiple stable steady states across temperatures including the transition 99 point. We generated 10,000 random initial conditions by scaling each chemical concentration of the 100 original steady state by 10௎ (ିଷ .଴ ,ଷ .଴ ). Here 𝑈(𝑎, 𝑏) denotes a uniform distribution ranging from 𝑎 101 to 𝑏. We sampled initial conditions uniformly in the logarithmic scale over the interval (10ିଷ , 10ଷ ). 102 See Materials and Methods for details of sampling procedure. Starting from each of these points, we 103 numerically integrated the system dynamics to 𝑡 = 10ଵ଴ 𝑠 and checked convergence for all cases. 104 This procedure revealed three coexisting stable steady states above 14.2 °C and two below 105 it. These results indicate that the discontinuous shift arises from a bifurcation, in which the 106 physiologically relevant steady state above 14.2 °C becomes destabilized as temperature decreases. 107 Bifurcation diagrams for all small molecules are provided in Supplementary Figure S1. Notably, the 108 disappearing steady state corresponds to the physiological E. coli state, characterized by aerobic 109 respiration and A TP production from glucose. 110 Examination of the steady state after the bifurcation revealed substantial changes in 111 metabolic fluxes and ATP production efficiency. The discontinuous shift caused glycolytic and TCA 112 cycle fluxes to decrease by factors of 100–1,000, while flux through the pentose phosphate pathway 113 (PPP) increased (Figure 2D). As a result, ATP yield—defined as the sum of ATP-producing fluxes 114 divided by the glucose uptake rate—declined following the temperature-induced bifurcation (Figure 115 2E). 116 To investigate the cause of the decreased ATP yield, we examined metabolite concentrations 117 along glycolysis and the TCA cycle. While most metabolite levels decreased after the temperature-118 induced bifurcation, the concentration of D-glycerate 1,3-diphosphate (13dpg) increased by two orders 119 of magnitude compared with its pre-bifurcation level (Figure 2C; Supplementary Figure S1). Despite 120 this elevation, the reaction degrading 13dpg—catalyzed by phosphoglycerate kinase (PGK), which 121 converts 13dpg and ADP to 3-phosphoglycerate and ATP in glycolysis (Figure 1)—showed a reduced 122 rate. This reduction was caused by a ~40-fold increase in the ATP/ADP ratio, driven by elevated ATP 123 and decreased ADP concentrations (Figure 2F; Supplementary Figure S1). Because the PGK reaction 124 lies near the center of glycolysis, its slowdown created a bottleneck that diminished flux through both 125 glycolysis and the downstream TCA cycle. 126 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint 127 128 Figure 2. (A) Distribution of reaction rate constants in the system. (B) Time courses of relative metabolite concentrations with all rate constants fixed at their values for 0ௗ°C and initial concentrations set to the original steady state; each line represents an individual metabolite. (C) Steady- state concentration of D-glycerate 1,3- bisphosphate at each temperature. (D) Relative reaction rates at steady state for each temperature, normalized to the rates at the original steady state. (E, F) Steady- state ATP yield (E) and A TP/ADP ratio (F) at each temperature. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint Fixing the ATP/ADP ratio prevents the temperature-induced bifurcation 129 The above results suggested that a sudden rise in the ATP/ADP ratio can trigger the temperature-130 induced bifurcation. To investigate this relationship, we maintained the A TP/ADP ratio nearly constant 131 by introducing a rapid, reversible ATP–ADP exchange reaction into the system. Because the ATP/ADP 132 ratio at the steady state before the bifurcation was approximately 1:1, the forward and reverse rate 133 constants were set to 10⁶, far exceeding the median values of other reactions in the model. A detailed 134 procedure for fixing the ATP/ADP ratio during temperature decrease is provided in Materials and 135 Methods. 136 The modified model with the ATP–ADP exchange reaction maintained a single, stable 137 steady state under the 37 ௗ°C parameter conditions. The resulting concentrations of 13dpg, the 138 ATP/ADP ratio, and the ATP yield were similar to those of the original steady state (Figure ௗ3A–C). 139 With the exchange reaction, lowering the temperature from 37ௗ°C to 0ௗ°C no longer caused a sharp rise 140 in 13dpg concentration or discontinuous changes in pathway fluxes, as shown by the relative changes 141 in reaction rates (Figureௗ3D). Consistently, the A TP/ADP ratio remained near unity, and the ATP yield 142 showed no abrupt decline across the temperature range (Figureௗ3B,ௗC). To test for multiple attractors 143 at each temperature, we generated 10,000 random initial points by scaling each metabolite 144 concentration of the original steady state by 10௎ (ିଷ .଴ ,ଷ .଴ ) and simulated the dynamics until 145 convergence. All trajectories converged to a single steady state at every temperature. These results 146 indicate that, if the ATP/ADP ratio is fixed at approximately 1.0, the system does not undergo 147 bifurcation, even though individual reaction rate constants still vary with temperature. 148 149 Maintaining the ATP yield by adjusting enzyme abundances 150 The analysis in the previous section demonstrated that the ATP yield can be maintained at a high level 151 by fixing the ATP/ADP ratio using a rapid, reversible ATP–ADP exchange reaction. However, cells do 152 not naturally possess such a reaction, and precise control of the ATP/ADP ratio in this manner is 153 biologically implausible. Instead, cells adjust intracellular enzyme abundances to adapt to diverse 154 environmental perturbations, including temperature fluctuations, osmotic stress, oxidative conditions, 155 and nutrient limitations (Jozefczuk et al., 2010; Schmidt et al., 2016; Ye et al., 2012; Zhao et al., 2024). 156 We therefore investigated whether ATP yield could be maintained under decreasing temperatures by 157 modulating enzyme abundances. To this end, we sampled steady states with varied enzyme levels 158 under low-temperature conditions and identified enzyme abundances that supported a high A TP yield. 159 160 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint 161 162 Figure 3. Changes in steady states of the modified model with A TP/ADP ratio homeostasis. (A) Steady-state concentration of D-glycerate 1,3-bisphosphate at each temperature. (B, C) Steady- state ATP/ADP ratio (B) and A TP yield (C) at each temperature. (D) Relative changes in reaction rates at steady state for each temperature. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint A naïve approach of randomly varying enzyme levels and integrating the system dynamics 163 until convergence for each set would require computationally intensive simulations. To overcome this, 164 we developed a faster method to sample numerous low-temperature steady states. Because every 165 reaction in our model is decomposed into elementary steps, each reaction rate function depends on at 166 most one or two variables and follows mass-action kinetics. Consequently, the rate functions are at 167 most quadratic polynomials, with the second variable, if present, corresponding to the concentration 168 of a small-molecule substrate. By fixing the concentrations of small molecules near their values in the 169 original steady state, the rate functions become linear, transforming the system of differential equations 170 into a linear one. Steady states of this linear system can then be obtained simply by solving a linear 171 equation (equationௗ(7) in Materials and Methods). 172 The linearized system has infinitely many steady states. To sample one, we first selected a 173 random reference point of enzyme abundances 𝑥௥௘௙ , near those of the original steady state. Each 174 concentration in the original steady state was multiplied by a value drawn from a log-normal 175 distribution over the interval (1/3, 3.0). Broadening this range reduces the fraction of sampled states 176 that are stable, so we chose the widest interval that still yielded stable solutions with reasonable 177 frequency. We then solved a linear optimization problem to identify the steady state closest to 𝑥௥௘௙ , 178 subject to additional constraints: (1) the glucose transport reaction must proceed inward while other 179 transport reactions act outward, and (2) the total enzyme amount must equal that of the original steady 180 state. These transport constraints were imposed to ensure the biological plausibility of the sampled 181 steady states. A mathematical formulation and detailed description of the sampling procedure are 182 provided in Materials and Methods. This approach samples steady states more efficiently than the 183 naïve method. Because stability was not enforced during sampling, we assessed it post hoc using the 184 Jacobian matrix at each steady state. Materials and Methods also describe how steady states were 185 classified as stable based on the Jacobian. 186 Using this method, we sampled 5,000 steady states with rate constants corresponding to 0 ௗ°C and 187 visualized the 94-dimensional enzyme profiles in two dimensions via principal component analysis 188 (PCA; Figureௗ4A). Among all computed solutions, the stable steady states were concentrated in a 189 specific region, where TCA cycle enzyme levels were elevated. Figure ௗ4B shows the sum of TCA 190 enzyme abundances at each steady state in the PCA plot. These results suggest that increasing TCA 191 cycle enzyme concentrations is important for stabilizing the metabolic state under low-temperature 192 conditions. 193 194 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint 195 196 Figure 4. (A) PCA plot of 94- dimensional enzyme abundances at steady states obtained with parameters for 0 ௗ°C. Explained variance: PC1, 16.9%; PC2, 4.7%. Stability of each steady state is indicated by color: blue, stable; grey, unstable. (B) TCA cycle enzyme abundances (color- coded) at steady states under low-temperature conditions. (C) Box plot of ATP yield for stable steady states at each temperature. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint The above results demonstrate that adjusting enzyme abundances allows the system to 197 achieve stable steady states near the original steady state, even under low-temperature conditions. 198 However, ATP yields at these stable states are not necessarily high. To evaluate A TP yield under low 199 temperatures, we sampled stable steady states across temperatures from 37ௗ°C to 0ௗ°C using the same 200

and assessed their A TP production. The A TP yield remained high even at low temperatures. 201 Figureௗ4C shows box plots of ATP yields for stable steady states at 37 ௗ°C, 19ௗ°C, and 0 ௗ°C. These 202 findings indicate that cells can maintain a high ATP yield during a low-temperature shift by modulating 203 enzyme abundances. 204 In practice, cells may not experience a sudden drop in temperature but rather undergo 205 gradual cooling. We therefore investigated which enzyme reallocations would allow the system to 206 maintain a high ATP yield as temperature decreases from 37 °C to 0ௗ°C. To do this, we fixed the target 207 concentrations of free enzymes and substrate-enzyme complexes ( 𝑥௥௘௙ ), and computed the steady-208 state enzyme abundances at each temperature. Along this path, enzyme levels at stable steady states 209 changed with temperature: TCA cycle enzymes increased, whereas ETC enzyme levels decreased. 210 Figuresௗ5A and 5B show the summed abundances of TCA cycle and ETC enzymes at stable steady 211 states across temperatures. This trend is consistent with experimental observations that cold-adapted 212 cells upregulate TCA cycle enzymes and downregulate ATP synthase subunits (Knapp et ௗal., 2024; 213 Phadtare & Inouye, 2004). 214 Because rate constants also vary with temperature, changes in enzyme abundance do not 215 translate directly into proportional changes in flux. Therefore, to understand how enzyme adjustments 216 maintain a high ATP yield during cooling, it is essential to examine the resulting metabolic fluxes 217 rather than enzyme levels alone. We analyzed reaction rates normalized by glucose influx to 218 investigate flux-level modulations arising from enzyme redistribution. The analysis revealed that, as 219 temperature decreased, normalized flux through the TCA cycle increased, along with fluxes of 220 reactions involved in ATP synthesis or consumption. Figuresௗ5C and 5D show the summed reaction 221 fluxes of the TCA cycle and ATP-utilizing reactions at stable steady states for each temperature. These 222 findings indicate that the adjusted enzyme abundances enhanced ATP–ADP turnover, analogous to the 223 rapid ATP–ADP exchange artificially imposed in the previous section. Together, these results suggest 224 that maintaining a high A TP/ADP ratio through increased cofactor turnover allows cells to sustain high 225 ATP yield even at low temperatures. 226 227 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint 228 Figure 5. (A, B) Changes in the summed abundances of TCA cycle enzymes (A) and ETC enzymes (B) at stable steady states across temperatures. Each line represents one 𝒙𝒓𝒆𝒇. (C, D) Changes in the summed absolute values of TCA cycle reaction rates (C) and A TP-utilizing reaction rates (D), both normalized to glucose influx, at stable steady states across temperatures. Each line represents one 𝒙𝒓𝒆𝒇. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint

229 In this study, we computationally investigated how perturbations in environmental parameters 230 influence the stability of a metabolic state. By varying rate constants in a temperature-dependent 231 manner, we found that lowering the temperature caused a bifurcation and a steep drop in the ATP yield. 232 This finding provides a model-based illustration of how complex metabolic systems can be 233 destabilized by parameter perturbations (Chakrabarti et al., 2013; Murayama et al., 2017; Steuer et al., 234 2006). 235 One way to maintain a high ATP yield was to keep the ATP/ADP ratio constant through 236 external control. Several studies have reported that coenzymes such as ATP and ADP exert a major 237 influence on the perturbation responses of kinetic models (Hatakeyama & Furusawa, 2017; Himeoka 238 & Mitarai 2022 ; Himeoka & Furusawa, 2025). In our model, fixing this cofactor ratio preserved a 239 stable metabolic state across the temperature range and maintained a high A TP yield as shown in Figure 240 3C. We also demonstrated that the A TP yield can be maintained at a high level by adjusting enzyme 241 abundances. Under such enzyme-tuning scenarios, the relative flux through A TP-utilizing reactions 242 increased, keeping the A TP/ADP ratio nearly constant. Thus, regulation of enzyme abundances can, in 243 principle, achieve the same effect as cofactor-ratio control. 244 Analysis of enzyme abundances to maintain the ATP yield revealed an increase in the TCA 245 cycle enzymes and a decrease in ETC enzymes. This trend aligns with experimental reports that E. 246 coli increases the TCA cycle enzyme levels and decreases ATP synthase subunits under cold 247 temperature (Knapp et al., 2024; Phadtare & Inouye, 2004) . This suggests that cells may modulate 248 enzyme abundances in vivo to maintain the ATP/ADP ratio and prevent an abrupt shift in their 249 metabolic state. Nevertheless, these findings require experimental validation, such as comprehensive 250 fluxomic analyses that cover pathways beyond our model. 251 Our analysis perturbed only temperature, starting from parameters fitted to the 37ௗ°C metabolic state. 252 Therefore, the conclusions depend on that particular parameter set. Additionally, the model is limited 253 to central metabolism and omits lipid pathways and many amino-acid biosynthesis routes that are 254 reported to be important for cold adaptation (Phadtare & Inouye, 2008; Zhang & Gross, 2021) . 255 Network structure modifications in kinetic models are known to exert strong influences on responses 256 for parameter perturbations (Hishida et al., 2023). Future work can probe cellular responses to cooling 257 with higher quantitative accuracy using a more detailed network and more precise temperature 258 dependence of parameters. 259 The same computational framework can be applied to environmental perturbations beyond 260 temperature. By allowing kinetic parameters to vary with factors such as pH, osmotic pressure, or 261 nutrient availability, we can reveal which metabolic functions become unstable and what operations—262 such as cofactor balance control or targeted enzyme tuning—can restore stability. 263 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint Overall, our results highlight the latent risk of destabilization in complex metabolic systems 264 under environmental perturbations and identify cofactor balance—specifically the ATP/ADP ratio—265 as a critical variable for maintaining stability. Applying this approach to other external parameters may 266 allow systematic identification of key variables that cells must regulate to preserve metabolic stability. 267 Such insights would advance our understanding of the principles underlying biological robustness. 268 269 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint 270 Supplementary Figure S1. Bifurcation diagram of small chemicals. The blue line represents the steady state used for parameter fitting in Khodayari etௗal. (2014) at 37ௗ°C, while the red and grey lines correspond to steady states that exist below 14.2 ௗ°C. The steady state at the blue line shifts to that represented by the red line due to the temperature-induced bifurcation. .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint 271 Reaction Pathway Formula 2DDA7Ptpp_ex 1 2dda7p 3PGtpp_ex 1 3pg ACALD 1 acald + 1 coa + 1 nad 1 accoa + 1 nadh ACALDtpp_ex 1 acald => ACK 1 ac + 1 atp 1 actp + 1 adp ACONTa TCA cycle 1 cit 1 acon-C + 1 h2o ACONTb TCA cycle 1 acon-C + 1 h2o 1 icit ACS 1 ac + 1 atp + 1 coa 1 accoa + 1 amp ACt2rpp_ex_H 1 ac => ADK 1 amp + 1 atp 2 adp AKGDH TCA cycle 1 akg + 1 coa + 1 nad 1 nadh + 1 succoa AKGt2rpp_ex 1 akg ALCD 1 etoh + 1 nad 1 acald + 1 nadh ATPM 1 atp 1 adp ATPS Electron transport chain 1 adp 1 atp CS TCA cycle 1 accoa + 1 oaa 1 cit + 1 coa CYTBD Electron transport chain 1 q8h2 1 q8 DDPA Pentose phosphate pathway 1 e4p + 1 pep 1 2dda7p DHAPtpp_ex 1 dhap D-LACt2pp_ex_H 1 lac-D => E4Ptpp_ex 1 e4p EDA 1 2ddg6p 1 g3p + 1 pyr .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint EDD 1 6pgc 1 2ddg6p ENO Glycolysis 1 2pg 1 pep ETOHt2rpp_ex_H 1 etoh => F6Ptpp_ex 1 f6p FBA Glycolysis 1 fdp 1 dhap + 1 g3p FBP Glycolysis 1 fdp 1 f6p FBP_Glpx 1 fdp 1 f6p FORt2pp_ex_H 1 for => FORtppi_ex 1 for => FRD TCA cycle 1 fum + 1 q8h2 1 q8 + 1 succ FRUpts2pp_ex 1 1 f6p + 1 pyr => 1 pep FUM TCA cycle 1 fum 1 mal FUMt2_2pp_ex_H 1 fum => G3Ptpp_ex 1 g3p G6PDH Pentose phosphate pathway 1 g6p + 1 nadp 1 6pgl + 1 nadph GAPD Glycolysis 1 g3p + 1 nad 1 13dpg + 1 nadh GLCptspp_ex_exi Glycolysis 1 pep 1 g6p + 1 pyr GLCS 1 adpglc 1 adp+ 1 glycogen GLGC 1 atp + 1 g1p 1 adpglc GLNS 1 atp + 1 glu 1 adp + 1 gln GLUD 1 glu + 1 nadp 1 akg + 1 nadph GLUN 1 gln 1 glu GLUS 1 akg + 1 gln + 1 nadph 2 glu + 1 nadp GLUt2rpp_ex_H 1 glu => .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint GND Pentose phosphate pathway 1 6pgc + 1 nadp 1 nadph + 1 ru5p HEX 1 atp + 1 glc-D 1 adp + 1 g6p ICDH TCA cycle 1 icit + 1 nadp 1 akg + 1 nadph ICL 1 icit 1 glx + 1 succ LDH 1 lac-D + 1 nad 1 nadh + 1 pyr MALS 1 accoa + 1 glx 1 coa + 1 mal MALt2_2pp_ex_H 1 mal => MDH TCA cycle 1 mal + 1 nad 1 nadh + 1 oaa MDH2 TCA cycle 1 mal + q8 oaa + q8h2 ME1 1 mal + 1 nad 1 nadh + 1 pyr ME2 1 mal + 1 nadp 1 nadph + 1 pyr NADH Electron transport chain 1 nadh + 1 q8 1 nad + 1 q8h2 NADTRHD 1 nad + 1 nadph 1 nadh + 1 nadp OAAtpp_ex 1 oaa PDH Glycolysis 1 coa + 1 nad + 1 pyr 1 accoa + 1 nadh PEPtpp_ex 1 pep PFKA Glycolysis 1 atp + 1 f6p 1 adp + 1 fdp PFKB Glycolysis 1 atp + 1 f6p 1 adp + 1 fdp PFL 1 coa + 1 pyr 1 accoa + 1 for PGI Glycolysis 1 g6p 1 f6p PGK Glycolysis 1 13dpg + 1 adp 1 3pg + 1 atp PGL Pentose phosphate pathway 1 6pgl 1 6pgc .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint PGL_spon Pentose phosphate pathway 1 6pgl 1 6pgc PGM Glycolysis 1 2pg 1 3pg PGMT 1 g1p 1 g6p PPC 1 pep 1 oaa PPCK 1 atp + 1 oaa 1 adp + 1 pep PPS Glycolysis 1 atp + 1 pyr 1 amp + 1 pep PTA 1 accoa 1 actp + 1 coa PYKA Glycolysis 1 adp + 1 pep 1 atp + 1 pyr PYKF Glycolysis 1 adp + 1 pep 1 atp + 1 pyr PYRpp_ex 1 pyr PYRt2rpp_ex_H 1 pyr R5PP Pentose phosphate pathway 1 r5p 1 rib RIB_Dtpp_ex 1 rib RPE Pentose phosphate pathway 1 ru5p 1 xu5p RPI Pentose phosphate pathway 1 r5p 1 ru5p RU5Ptpp_ex 1 ru5p => S7Ptpp_ex 1 s7p => SUCCt2_2pp_ex_H 1 succ => SUCCt3pp_ex_H 1 succ SUCD TCA cycle 1 q8 + 1 succ 1 fum + 1 q8h2 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint SUCOAS TCA cycle 1 atp + 1 coa + 1 succ 1 adp + 1 succoa TALA Pentose phosphate pathway 1 g3p + 1 s7p 1 e4p + 1 f6p THD 1 nadh + 1 nadp 1 nad + 1 nadph TKT1 Pentose phosphate pathway 1 r5p + 1 xu5p 1 g3p + 1 s7p TKT2 Pentose phosphate pathway 1 e4p + 1 xu5p 1 f6p + 1 g3p TPI Glycolysis 1 dhap 1 g3p 272 Supplementary Table S1. List of all reactions in the model used in our simulations. 273 274 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint Reaction Regulator Regulation Type ACALD amp Inhibition EDA 6pgc Inhibition EDA g3p Inhibition FBA dhap Inhibition FBP adp Inhibition FBP amp Inhibition FBP pi Inhibition FBP pep Inhibition FUM cit Inhibition G6PDH nadh Inhibition G6PDH nadph Inhibition GLGC amp Inhibition GLUN adp Inhibition GLUN atp Inhibition GLUN glu Inhibition GLUN nad Inhibition GLUN nadh Inhibition GLUN nadp Inhibition GLUN nadph Inhibition GLUN pi Inhibition GLUS glu Inhibition GND atp Inhibition GND ru5p Inhibition GND fdp Inhibition GND nadp Inhibition ICDH glx Inhibition ICDH oxa Inhibition .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint ICDH pep Inhibition ICL akg Inhibition ICL acon Inhibition ICL pep Inhibition ICL succ Inhibition LDH atp Inhibition PGI 6pgc Inhibition PGI pep Inhibition PGMT accoa Inhibition PGMT g1p Inhibition PPC mal Inhibition PPC pep Inhibition PPCK atp Inhibition PPCK pep Inhibition PPS akg Inhibition PPS adp Inhibition PPS amp Inhibition PPS oaa Inhibition PPS pep Inhibition PTA adp Inhibition PTA atp Inhibition PTA nadh Inhibition PTA nadph Inhibition PTA pep Inhibition PTA pyr Inhibition PYKA PYKF atp Inhibition PYKA PYKF succoa Inhibition RPI amp Inhibition .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint TALA pi Inhibition 275 Supplementary Table S2. List of all regulations in the model used in our simulations. 276 277 278 279 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint

280 281 Model Description 282 The model used in this analysis is based on the model from Himeoka et al. (2025), which is a modified 283 version of the Khodayari et al. (2014) model, adjusted to ensure the existence of a stable steady state. 284 In the present study, the parameters of inflow and outflow reactions were further modified so that 285 glucose serves as the sole carbon source in the temperature range used for the simulations. The model 286 includes glycolysis, the TCA cycle, the pentose phosphate pathway (PPP), glutamine synthesis, the 287 electron transport chain (ETC), and the Entner-Doudoroff pathway. See Supplementary Table S1 and 288 S2 for complete reaction formulae, pathway components, and regulations. 289 In our model, all reactions are decomposed into elementary reactions that represent 290 substrate-enzyme binding and dissociation. This is illustrated using a simple reversible reaction as 291 follows: 292 For a reversible reaction A + B ⇌ C catalyzed by an enzyme E, we write: 293 A + E ⇌ EA 294 EA + B ⇌ EAB 295 EAB ⇌ EC 296 EC ⇌ E + C. 297 All steps are reversible, and the forward and reverse directions are treated as distinct elementary 298 reactions. The total concentration of free and complex enzymes does not change in time: 299 [𝐸] + [𝐸𝐴] + [𝐸𝐴𝐵] + [𝐸𝐶] = const. 303 The rates of these elementary reactions are calculated based on mass-action kinetics. Note that in our 300 model, chemical concentrations are dimensionless because they are normalized to 1.0 at the reference 301 steady state as a standard procedure in ensemble modeling (Khodayari et al., 2014; Tran et al., 2008). 302 In our model, the time evolution of concentrations is described by an ODE system 304 𝑑𝑥(𝑡) 𝑑𝑡 = 𝑆𝑟൫𝑥(𝑡); 𝑘(𝑇)൯, (1) 311 where 𝑥 is the concentration vector of small molecules, free enzymes, and substrate-enzyme 305 complexes; 𝑆 is the stoichiometric matrix; 𝑘(𝑇) is the set of rate constants at temperature 𝑇; and 306 𝑟(𝑥) is the vector of reaction rate functions. All reactions are decomposed into elementary reactions, 307 and each reaction rate is modeled by the law of mass action 308 𝑟௜(𝑥) = 𝑘௜(𝑇) ෑ 𝑥௝ ఔ ೕ೔ ௡ ௝ୀଵ , (2) 312 with exponents 𝜈௝௜ ∈ {0,1} determined by the stoichiometry and 𝑘௜(𝑇) the temperature-dependent 309 rate constants (specified in the next section). 310 ODE simulations were performed using the ode15s function in MATLAB (The MathWorks 313 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint Inc., 2023). 314 315 Temperature-Dependent Rate Constants 316 Rate constants were adjusted with the Eyring-Polanyi equation 317 𝑘(𝑇) = 𝑘஻ 𝑇 ℎ exp ቆ − ∆𝐺‡ 𝑅𝑇 ቇ (3) 324 for each temperature. Here, 𝑘 is the rate constant, ℎ is the Planck constant, 𝑘஻ is the Boltzmann 318 constant, 𝑅 is the gas constant, and 𝑇 is the absolute temperature (𝐾). The activation free energy 319 ∆𝐺‡ for each reaction is back-calculated from parameter values obtained in Khodayari et al. 2014 320 using the following equation: 321 ∆𝐺‡ = −𝑅𝑇௥௘௙ ln ቆ 𝑘∗ ℎ 𝑘஻ 𝑇௥௘௙ ቇ , (4) 325 with 𝑇௥௘௙ = 310.15, and 𝑘∗ is the rate constant at 𝑇௥௘௙ . Based on this value, reaction rate constants 322 at altered temperatures were calculated. 323 326 Maintaining the ATP/ADP Ratio 327 To forcibly maintain the ATP/ADP ratio while lowering the temperature, we added a new reversible 328 reaction that exchanges A TP and ADP to the model. The rate constant for this exchange was set to 329 10଺ , a value larger than the median of all rate constants at 37 °C. During low-temperature simulations, 330 this constant was not adjusted for temperature dependence. The concentrations of the free metabolites 331 at the steady state are normalized to unity in the ensemble modeling procedure carried out for the 332 model construction in Khodayari et al. (2014). Thus, ATP/ADP ratio in the original model is 1:1 in the 333 scale of normalized concentrations. Therefore, the force-balancing A TP/ADP ratio to unity 334 corresponds to the recovery of the ATP/ADP ratio in the original model. 335 336 Perturbed Initial Condition Sampling Procedure 337 To identify multiple stable steady states in the system, we perform numerical integrations starting from 338 perturbed initial conditions with specified temperature parameters. These initial conditions are 339 generated by multiplying the concentrations of an original steady state by random values with ensuring 340 that the conserved quantities are not changed. 341 If a vector 𝑑 satisfies 342 𝑑୘ 𝑆 = 0, 345 then the inner product of 𝑑 and a concentration vector 𝑥 does not change in time, since 343 𝑑 𝑑𝑡 (𝑑୘ 𝑥) = 0. 346 344 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint If random perturbations are applied directly to the concentrations of the original steady state 347 𝑥௢௥௚ , the resulting concentration vector 𝑥௥௔௡ௗ will generally have conserved quantities different from 348 those of 𝑥௢௥௚ . To obtain perturbed concentrations that preserve the conserved quantities of 𝑥௢௥௚ , we 349 search for the point closest to 𝑥௥௔௡ௗ that satisfies the conservation constraints. 350 This is formulated as the following optimization problem: 351 𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 ‖𝑥 − 𝑥௥௔௡ௗ ‖ଵ 354 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 355 𝑑௜ ୘ 𝑥 = 𝑑௜ ୘ 𝑥௢௥௚ , for all 𝑖. 356 Where {𝑑௜} form a basis of ker 𝑆୘ . This optimization is done using Gurobi Optimizer (Gurobi 352 Optimization, 2024). 353 357 Steady State Sampling Procedure 358 To identify enzyme abundances that yield a stable metabolic state with low-temperature parameters, 359 we sampled steady states for each temperature-specific parameter set and computed the corresponding 360 enzyme abundances. Each reaction is decomposed into elementary reactions whose reactants are (i) 361 substrate and enzyme, (ii) substrate and substrate-enzyme complex, or (iii) substrate-enzyme complex 362 only. Therefore, each reaction rate function is a first- or second-order monomial in its variables. 363 Fixing concentrations of substrates (small molecules) at their values near the original steady 364 state turns every rate into a linear function of the free-enzyme or substrate-enzyme complex 365 concentrations. 𝑥௘௡௭௬௠௘ is a vector of the remaining variables, i.e., the concentrations of free 366 enzymes and substrate-enzyme complexes. Using a matrix 𝑉, whose non-zero values are determined 367 by the temperature-specific rate constants and the fixed substrate concentrations, we have 368 𝑟 = 𝑉𝑥௘௡௭௬௠௘ (5) 370 where 𝑟 is the vector of reaction rates. 369 For example, if the system has only one reaction A + B ⇌ C catalyzed by an enzyme E, 371 reaction rate functions are 372 𝑟 = ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ 𝑘ଵ௙ [𝐴][𝐸] 𝑘ଵ௕ [𝐸𝐴] 𝑘ଶ௙ [𝐸𝐴][𝐵] 𝑘ଶ௕ [𝐸𝐴𝐵] 𝑘ଷ௙ [𝐸𝐴𝐵] 𝑘ଷ௕ [𝐸𝐶] 𝑘ସ௙ [𝐸𝐶] 𝑘ସ௕ [𝐸][𝐶] ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞ , 375 where 𝑘௜ is rate constant of forward or backward direction of each elementary reaction. When the 373 substrate concentrations are fixed at ([𝐴]∗, [𝐵]∗, [𝐶]∗)୘ , the rate functions can be rewritten as 374 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint 𝑟 = ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ 𝑘ଵ௙ [𝐴]∗ 0 0 0 0 𝑘ଵ௕ 0 0 0 𝑘ଶ௙ [𝐵]∗ 0 0 0 0 𝑘ଶ௕ 0 0 0 𝑘ଷ௙ 0 0 0 0 𝑘ଷ௕ 0 0 0 𝑘ସ௙ 𝑘ସ௕ [𝐶]∗ 0 0 0 ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞ ൮ [𝐸] [𝐸𝐴] [𝐸𝐴𝐵] [𝐸𝐶] ൲ . 376 Using submatrix of the stoichiometric matrix corresponding to the remaining variables, 377 denoted as 𝑆௘௡௭௬௠௘ , the time-course change of 𝑥௘௡௭௬௠௘ can be written as 378 𝑑𝑥௘௡௭௬௠௘ 𝑑𝑡 = 𝑆௘௡௭௬௠௘ 𝑉𝑥௘௡௭௬௠௘ . (6) 380 Therefore, steady-state solutions are obtained just by solving 379 𝑆௘௡௭௬௠௘ 𝑉𝑥௘௡௭௬௠௘ = 0. (7) 381 Eq. (7) is less-determined linear equation when 𝑆௘௡௭௬௠௘ 𝑉 is not full-rank. In that case, Eq. 382 (7) has infinitely many solutions. Among such solutions of Eq. (7), we select the one closest to a 383

point ൫𝑥௥௘௙ ൯, which is randomly sampled around the enzyme abundances of the original 384 steady state. We also require steady states that 385 1. Glucose is a single energy resource, and transport reactions of other chemicals work as efflux. 386 2. The total enzyme abundance of all reactions is equal to that of the original steady state. 387 Steady states satisfying those conditions are obtained by solving the following optimization problem: 388 𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 ‖𝑥௘௡௭௬௠௘ − 𝑥௥௘௙ ‖ଵ 392 𝑆𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 393 𝑆௘௡௭௬௠௘ 𝑉𝑥௘௡௭௬௠௘ = 0 394 ෍ 𝑥௘௡௭௬௠௘ ,௜ ௜ = const. 395 𝑣ீ௟௨௖௢௦௘ூ௡௣௨௧ > 0 396 𝑣ை௧௛௘௥்௥௔௡௦௣ < 0 397 The solution 𝑥௘௡௭௬௠௘ obtained from this optimization, together with the fixed small molecule 389 concentrations near the 37 °C steady state, yields a steady state of the full system. These optimizations 390 are done using Gurobi Optimizer (Gurobi Optimization, 2024). 391 398 Stability Assessment of Steady States 399 The stability of each steady state found by the optimization was determined by computing the 400 eigenvalues of the system’s Jacobi matrix at the steady state. The model contains 98 conserved 401 quantities, so the full 643 × 643 Jacobi matrix has 98 zero eigenvalues. The zero eigenvalues are often 402 tricky to distinguish from the very small, non-zero eigenvalues due to numerical errors. Thus, we 403 avoided obtaining zero eigenvalues by pre-solving the dependent variables using the conserved 404 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint quantities. Note that one conserved quantity leads to a single linear equation. For each conserved 405 quantity, we choose a single chemical species and solve for its concentration as a function of the 406 remaining species’ concentrations. From this pre-solved ODE system, we obtain the reduced 545 × 407 545 Jacobi matrix. A steady state was classified as unstable if the largest real part of any eigenvalue 408 of the reduced Jacobi matrix was positive. All eigenvalues were computed with interval-arithmetic 409 routines in INTLAB, ensuring that numerical errors did not flip the sign of any eigenvalue (Rump, 410 1999). 411 412

413 This work is supported in part by JSPS KAKENHI (23KJ1324 to A . H.; 22K15069, 24H01118 and 414 25H01390 to Y . H.; 22K21344 and 23H0247 to C. F.), and GteX Program Japan Grant No. 415 JPMJGX23B4. 416 417 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint

418 Arcus, V . L., & Mulholland, A. J. (2020). Temperature, Dynamics, and Enzyme-Catalyzed Reaction 419 Rates. Annual Review of Biophysics, 49(1), 163–180. https://doi.org/10.1146/annurev-biophys-420 121219-081520 421 Chakrabarti, A., Miskovic, L., Soh, K. C., & Hatzimanikatis, V . (2013). Towards kinetic modeling of 422 genome-scale metabolic networks without sacrificing stoichiometric, thermodynamic and 423 physiological constraints. Biotechnology Journal , 8(9), 1043–1057. 424 https://doi.org/10.1002/biot.201300091 425 Christodoulou, D., Link, H., Fuhrer, T., Kochanowski, K., Gerosa, L., & Sauer, U. (2018). Reserve 426 Flux Capacity in the Pentose Phosphate Pathway Enables Escherichia coli’s Rapid Response to 427 Oxidative Stress. Cell Systems, 6(5), 569-578.e7. https://doi.org/10.1016/j.cels.2018.04.009 428 El-Samad, H., Kurata, H., Doyle, J. C., Gross, C. A., & Khammash, M. (2005). Surviving heat shock: 429 Control strategies for robustness and performance. Proceedings of the National Academy of 430 Sciences, 102(8), 2736–2741. https://doi.org/10.1073/pnas.0403510102 431 Evans, M. G., & Polanyi, M. (1935). Some applications of the transition state method to the calculation 432 of reaction velocities, especially in solution. Transactions of the Faraday Society , 31, 875. 433 https://doi.org/10.1039/tf9353100875 434 Eyring, H. (1935). The Activated Complex in Chemical Reactions. The Journal of Chemical Physics, 435 3(2), 107–115. https://doi.org/10.1063/1.1749604 436 Gruber, C. H., Noor, E., Buffing, M. F., & Sauer, U. (2025). Systematic identification of allosteric 437 effectors in Escherichia coli metabolism. Proceedings of the National Academy of Sciences of 438 the United States of America, 122(10). https://doi.org/10.1073/pnas.2423767122 439 Gurobi Optimization, L. (2024). Gurobi Optimizer Reference Manual. 440 Hatakeyama, T. S., & Furusawa, C. (2017). Metabolic dynamics restricted by conserved carriers: 441 Jamming and feedback. PLoS Computational Biology , 13(11). 442 https://doi.org/10.1371/journal.pcbi.1005847 443 Himeoka, Y ., & Mitarai, N. (2022). Emergence of Growth and Dormancy from a Kinetic Model of the 444 Escherichia Coli Central Carbon Metabolism., Physical Review Research 4 (4): 043223. 445 Himeoka, Y ., & Furusawa, C. (2025). Perturbation-response analysis of in silico metabolic dynamics 446 revealed hard-coded responsiveness in the cofactors and network sparsity. ELife, 13. 447 https://doi.org/10.7554/eLife.98800 448 Hishida, A., Okada, T., & Mochizuki, A. (2023). Patterns of change in regulatory modules of chemical 449 reaction systems induced by network modification. PNAS Nexus , 3(1). 450 https://doi.org/10.1093/pnasnexus/pgad441 451 Jozefczuk, S., Klie, S., Catchpole, G., Szymanski, J., Cuadros-Inostroza, A., Steinhauser, D., Selbig, 452 J., & Willmitzer, L. (2010). Metabolomic and transcriptomic stress response of Escherichia coli. 453 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint Molecular Systems Biology, 6. https://doi.org/10.1038/msb.2010.18 454 Khodayari, A., Zomorrodi, A. R., Liao, J. C., & Maranas, C. D. (2014). A kinetic model of Escherichia 455 coli core metabolism satisfying multiple sets of mutant flux data. Metabolic Engineering, 25, 456 50–62. https://doi.org/10.1016/j.ymben.2014.05.014 457 Knapp, B. D., Willis, L., Gonzalez, C., Vashistha, H., Jammal-Touma, J., Tikhonov, M., Ram, J., 458 Salman, H., Elias, J. E., & Huang, K. C. (2024). Metabolic rearrangement enables adaptation of 459 microbial growth rate to temperature shifts. Nature Microbiology , 10(1), 185–201. 460 https://doi.org/10.1038/s41564-024-01841-4 461 Miyakoshi, M. (2024). Multilayered regulation of amino acid metabolism in Escherichia coli. Current 462 Opinion in Microbiology, 77, 102406. https://doi.org/10.1016/j.mib.2023.102406 463 Murayama, Y ., Kori, H., Oshima, C., Kondo, T., Iwasaki, H., & Ito, H. (2017). Low temperature 464 nullifies the circadian clock in cyanobacteria through Hopf bifurcation. Proceedings of the 465 National Academy of Sciences of the United States of America , 114(22), 5641–5646. 466 https://doi.org/10.1073/pnas.1620378114 467 Noll, P., Lilge, L., Hausmann, R., & Henkel, M. (2020). Modeling and Exploiting Microbial 468 Temperature Response. Processes, 8(1), 121. https://doi.org/10.3390/pr8010121 469 Phadtare, S. (2004). Recent Developments in Bacterial Cold-Shock Response. Current Issues in 470 Molecular Biology. https://doi.org/10.21775/cimb.006.125 471 Phadtare, S., & Inouye, M. (2004). Genome-wide transcriptional analysis of the cold shock response 472 in wild-type and cold-sensitive, quadruple-csp-deletion strains of Escherichia coli. Journal of 473 Bacteriology, 186(20), 7007–7014. https://doi.org/10.1128/JB.186.20.7007-7014.2004 474 Phadtare, S., & Inouye, M. (2008). The Cold Shock Response. EcoSal Plus , 3(1). 475 https://doi.org/10.1128/ecosalplus.5.4.2 476 Riccardi, C., Calvanese, M., Ghini, V ., Alonso-Vásquez, T., Perrin, E., Turano, P., Giurato, G., Weisz, 477 A., Parrilli, E., Tutino, M. L., & Fondi, M. (2023). Metabolic Robustness to Growth Temperature 478 of a Cold- Adapted Marine Bacterium. MSystems, 8(2). 479 https://doi.org/10.1128/msystems.01124-22 480 Rump, S. M. (1999). INTLAB — INTerval LABoratory. In Developments in Reliable Computing (pp. 481 77–104). Springer Netherlands. https://doi.org/10.1007/978-94-017-1247-7_7 482 Sawers, R. G. (2025). How FocA facilitates fermentation and respiration of formate by Escherichia 483 coli. In Journal of Bacteriology (V ol. 207, Issue 2). American Society for Microbiology. 484 https://doi.org/10.1128/jb.00502-24 485 Schmidt, A., Kochanowski, K., Vedelaar, S., Ahrné, E., V olkmer, B., Callipo, L., Knoops, K., Bauer, 486 M., Aebersold, R., & Heinemann, M. (2016). The quantitative and condition-dependent 487 Escherichia coli proteome. Nature Biotechnology , 34(1), 104–110. 488 https://doi.org/10.1038/nbt.3418 489 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint Stephanopoulos, G. ., Aristidou, A. A. ., & Nielsen, J. Høiriis. (1998). Metabolic engineeringࣟ :490 principles and methodologies. Academic Press. 491 Steuer, R., Gross, T., Selbig, J., & Blasius, B. (2006). Structural kinetic modeling of metabolic 492 networks. Proceedings of the National Academy of Sciences , 103(32), 11868–11873. 493 https://doi.org/10.1073/pnas.0600013103 494 Strutz, J., Martin, J., Greene, J., Broadbelt, L., & Tyo, K. (2019). Metabolic kinetic modeling provides 495 insight into complex biological questions, but hurdles remain. Current Opinion in Biotechnology, 496 59, 24–30. https://doi.org/10.1016/j.copbio.2019.02.005 497 Sudarsan, S., Demling, P., Ozdemir, E., Ben Ammar, A., Mennicken, P., Buescher, J. M., Meurer, G., 498 Ebert, B. E., & Blank, L. M. (2025). Acetol biosynthesis enables NADPH balance during 499 nitrogen limitation in engineered Escherichia coli. Microbial Cell Factories , 24(1). 500 https://doi.org/10.1186/s12934-025-02687-z 501 The MathWorks Inc. (2023). MATLAB version: R2023b. The MathWorks Inc. 502 Tran, L. M., Rizk, M. L., & Liao, J. C. (2008). Ensemble modeling of metabolic networks. Biophysical 503 Journal, 95(12), 5606–5617. https://doi.org/10.1529/biophysj.108.135442 504 Wendering, P., & Nikoloski, Z. (2023). Model-driven insights into the effects of temperature on 505 metabolism. Biotechnology Advances , 67, 108203. 506 https://doi.org/10.1016/j.biotechadv.2023.108203 507 Ye, Y ., Zhang, L., Hao, F., Zhang, J., Wang, Y ., & Tang, H. (2012). Global metabolomic responses of 508 escherichia coli to heat stress. Journal of Proteome Research , 11(4), 2559–2566. 509 https://doi.org/10.1021/pr3000128 510 Zhang, Y ., & Gross, C. A. (2021). Cold Shock Response in Bacteria. Annual Review of Genetics, 55(1), 511 377–400. https://doi.org/10.1146/annurev-genet-071819-031654 512 Zhao, J., Chen, K., Palsson, B. O., & Yang, L. (2024). StressME: Unified computing framework of 513 Escherichia coli metabolism, gene expression, and stress responses. PLoS Computational 514 Biology, 20(2 February). https://doi.org/10.1371/journal.pcbi.1011865 515 516 517 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted September 1, 2025. ; https://doi.org/10.1101/2025.08.27.672763doi: bioRxiv preprint