A Kinetic Platform to Determine the Fate of Hydrogen Peroxide in Escherichia coli

Hydrogen peroxide (H2O2) is used by phagocytic cells of the innate immune response to kill engulfed bacteria. H2O2 diffuses freely into bacteria, where it can wreak havoc on sensitive biomolecules if it is not rapidly detoxified. Accordingly, bacteria have evolved numerous systems to defend themselves against H2O2, and the importance of these systems to pathogenesis has been substantiated by the many bacteria that require them to establish or sustain infections. The kinetic competition for H2O2 within bacteria is complex, which suggests that quantitative models will improve interpretation and prediction of network behavior. To date, such models have been of limited scope, and this inspired us to construct a quantitative, systems-level model of H2O2 detoxification in Escherichia coli that includes detoxification enzymes, H2O2-dependent transcriptional regulation, enzyme degradation, the Fenton reaction and damage caused by •OH, oxidation of biomolecules by H2O2, and repair processes. After using an iterative computational and experimental procedure to train the model, we leveraged it to predict how H2O2 detoxification would change in response to an environmental perturbation that pathogens encounter within host phagosomes, carbon source deprivation, which leads to translational inhibition and limited availability of NADH. We found that the model accurately predicted that NADH depletion would delay clearance at low H2O2 concentrations and that detoxification at higher concentrations would resemble that of carbon-replete conditions. These results suggest that protein synthesis during bolus H2O2 stress does not affect clearance dynamics and that access to catabolites only matters at low H2O2 concentrations. We anticipate that this model will serve as a computational tool for the quantitative exploration and dissection of oxidative stress in bacteria, and that the model and methods used to develop it will provide important templates for the generation of comparable models for other bacterial species.

As a result of its importance as a signaling molecule, the most complete models of H 2 O 2 metabolism currently available were developed for mammalian systems [22][23][24]. They have included H 2 O 2 elimination by antioxidants (e.g., glutathione and thioredoxin) and enzymes (e.g., catalase, glutathione peroxidase, glutathione reductase, glutaredoxin, and peroxiredoxin), and processing of oxidized protein thiols [22][23][24]. However, these models were specific to mammalian physiology and did not include transcriptional regulation, enzyme degradation, side reactions of H 2 O 2 with sensitive biomolecules such as methionine and pyruvate, or the related reactive oxygen species O 2 − • and •OH and their associated reactivity (e.g., •OH rapidly oxidizes all twenty amino acids and glutathione). Although models equivalent to those of mammalian systems have yet to be described for bacteria, there has been progress in modeling subsystems of the H 2 O 2 response network under H 2 O 2 stress, such as the thioredoxin system in E. coli [25].
Here, we have generated a kinetic model of H 2 O 2 stress in E. coli whose components are depicted in Fig 1. The biochemical reaction network is compartmentalized into media and intracellular spaces, includes spontaneous and enzymatic detoxification of H 2 O 2 ,  Table) and enzymatic reactions (blue lines, details in S3 Table). Metabolite abbreviations can be found in S1 Table. Information regarding enzyme degradation can be found in S2 Table, and enzyme expression is described in S3 Table. For clarity, •OH reaction products with amino acids and protons are not included in the diagram. transcriptional regulation and inactivation of detoxification enzymes, and reactions of H 2 O 2 and its degradation intermediates (e.g., •OH) with biomolecules (e.g., pyruvate, glutathione, all twenty amino acids). Parameters were informed from literature or trained using an iterative and integrated computational and experimental approach (Fig 2). The design criteria we chose to use to develop the model stipulated that consistent discrimination between clearance contributions by the major detoxification systems (AHP, HPI, and HPII) needed to be achieved. Once the design criteria were met, remaining parametric uncertainty was accounted for with use of a Markov chain Monte Carlo (MCMC) procedure to explore the viable parameter space and assemble an ensemble of models that performed comparably well with the training data [26]. The ensemble was then used to quantitatively investigate the importance of carbon availability and translation to H 2 O 2 detoxification, and its predictions were experimentally confirmed. Uncertain parameters in each of the ten model structures are optimized on wild-type clearance data of 10, 25, 100, and 400 μM H 2 O 2 , starting from 1,000 random initial parameter sets. Any models within an evidence ratio of 10 (ER 10) are used to calculate cumulative H 2 O 2 clearance by the different detoxification pathways. If the calculated H 2 O 2 distributions between the models are inconsistent, simulations are used to suggest experiments that differentiate between the disagreeing models, experiments are executed, and the optimization is performed on all experimental data for the model structures that had at least one parameter set with an ER 10. Once consistent H 2 O 2 distributions are realized, we identify an ensemble of parameter sets that can all describe the data comparably well with an MCMC procedure. We then assess whether H 2 O 2 distributions are consistent across the entire ensemble. If the calculated H 2 O 2 distributions are consistent, we undertake forward predictions. doi:10.1371/journal.pcbi.1004562.g002

Design criteria
Our aim was to construct a systems-level kinetic model of H 2 O 2 detoxification in E. coli that could provide consistent predictions of H 2 O 2 distributions among its different detoxification pathways after exposure to a range of initial H 2 O 2 boluses. To accomplish this goal in the most efficient way possible, we adopted the systematic approach shown in Fig 2. Briefly, we began with a minimal number of experiments, wild-type clearance of different initial H 2 O 2 concentrations. After optimizing uncertain parameters, we selected models based on their relative likelihood, also referred to as their evidence ratio (ER) [27][28][29][30][31][32], discarding models more than ten times less likely than the most-likely model in our set (ER!10). If the acceptable models did not uniformly attribute H 2 O 2 detoxification to the same pathways, we performed simulations to suggest experiments that could resolve the disagreement. Those experiments were then performed, and data used to arrive at updated parameter estimates. This process was continued until we arrived at a model or set of models that rendered consistent H 2 O 2 distributions. Since some parameters may not have been important to H 2 O 2 clearance under the conditions used here, and therefore, unlikely to be informed by the training procedure, we explored the parameter space using a previously developed MCMC procedure [26] to assemble an ensemble of parameter sets that could all describe the H 2 O 2 clearance data comparably well (ER 10). In this way, we could ensure that forward predictions were not dependent on ill-defined parameters. We note that this procedure also accounts for cases in which parameter pairings or more complex relationships rather than absolute values are important by varying all parameters simultaneously when walking away from known viable points. Also, before proceeding to forward predictions, we confirmed that all the models within the ensemble still satisfied the design criteria.

Kinetic model of H 2 O 2 metabolism in E. coli
We constructed a compartmentalized reaction network that includes spontaneous and enzymatic reactions present in an E. coli culture under H 2 O 2 stress, transcriptional regulation of AHP and HPI, and degradation/inactivation of the major detoxification enzymes AHP, HPI, and HPII. Uncertainty exists with regard to the dynamics of enzyme degradation/inactivation in the presence of H 2 O 2 , as well as the possibility of an H 2 O 2 gradient across the membrane. Specifically, the enzymes could be degraded or inactivated in an H 2 O 2 -independent manner, either with a fixed degradation constant [33,34], or optimized to account for the varying degradation rates of different proteins [35,36]. Alternatively, the H 2 O 2 detoxification enzymes could be poisoned by their own substrate [37][38][39], following bimolecular [40] or more complex kinetics [37].
In addition to the indeterminacy in degradation/inactivation kinetics, there is evidence supporting [41] and opposing [42] the presence of an H 2 O 2 gradient across the cell membrane. Models accounting for these various possibilities are presented in Table 1, along with their corresponding number of uncertain parameters. The introduction of unknown parameters has the potential to improve agreement between model simulations and experimental data solely by increasing the flexibility of the model. For this reason, we calculated the relative likelihood of models, otherwise known as their evidence ratio (ER) [27][28][29][30][31][32] based on their respective Akaike Information Criterion (AIC), which is a commonly used statistical metric that weighs goodness of fit against model complexity when discriminating between competing models [27,[43][44][45]. Models with a relative likelihood of ten times less than the best model in the set (ER 10) were considered acceptable, whereas others were discarded.

Experimentally-driven model refinement and parameter estimation
When optimizing parameters simultaneously on clearance of 10, 25, 100, and 400 μM boluses in wild-type cultures, 35 of the 10,000 models had an ER 10 and were considered viable models. The windows of simulation results of these 35 models are presented in Fig 3A-3D along with the experimental clearance data they were trained on. None of these models contained a gradient; 30 were structure 2 models, and 5 were structure 3. We note that structures that contain a gradient could be favored over those with no gradient under different conditions (e.g., training on data from a single H 2 O 2 concentration); however, our goal was to arrive at a model that could describe a wide range of bolus concentrations, and the gain in simulation accuracy for gradient models did not justify addition of the extra parameter as determined by the ER for the experimental conditions considered here. When the H 2 O 2 distributions of the acceptable models were analyzed, the utility of AHP and the catalases separated into two distinct groups at all bolus concentrations (Fig 3E-3H). The reaction fluxes through AHP and HPI+HPII can be found in S2A-S2D Fig, and those also separated into two distinct groups. We found that these two groups represented predictions made by the two model structures. At all bolus concentrations, structure 2 models predicted a greater contribution by AHP than did structure 3 models. Indeterminacy in catalase null mutant simulations (Fig 3I-3L) suggested that experiments on a strain lacking both HPI (katG) and HPII (katE) would resolve this discrepancy.
Simultaneous training of models on wild-type and ΔkatE ΔkatG clearance data was able to resolve the uncertainty between structures 2 and 3. This training iteration resulted in 965 models that all had an ER 10 ( Fig 4A-4D), and all acceptable models were structure 3, which suggests that bimolecular H 2 O 2 -dependent enzyme degradation is an important feature of the detoxification network. We note that clearance of 400 μM H 2 O 2 by ΔkatE ΔkatG was omitted because significant cell death was observed (S1H Fig), and the models were not designed to simulate cell death and possible lysis. All models predicted similar distributions across the major pathways (Fig 4E-4H), but diverged when we looked more closely at the individual catalase contributions (Fig 4I-4L). Reaction fluxes through the major pathways and individual catalases can be found in S2E- S2H Fig and S3A-S3D Fig, respectively. The different parameter sets predicted a range of clearance profiles after removal of either catalase (S4 Fig), suggesting that data obtained from these mutants would resolve the disagreement between models.
Training uncertain model parameters on wild-type, ΔkatE ΔkatG, ΔkatE, and ΔkatG data resulted in 40 parameters sets from the 1,000 random initializations that were within an ER of 10 ( Fig 5A-5D). All of these models agreed regarding how H 2 O 2 distributes across not only the major pathways (Fig 5E-5H), but also the individual catalases (Fig 5I-5L). Reaction fluxes  . Structure 2 (green) and structure 3 (purple) models predict different clearance dynamics after this perturbation, suggesting that data obtained from this mutant could be used to discriminate between the model structures.

Exploration of the viable parameter space
The identification of universal "sloppiness" in computational biological models [46], meaning many parameters are poorly constrained after fitting on experimental data, led to the development of a number of methods designed to identify ensembles of parameter sets that could comparably describe the data and be used to assess the robustness of forward predictions [26,46]. Methods such as "brute force" uniform sampling or Gaussian sampling become impossible with increasingly complex models, so computational biologists have turned to the use of Monte Carlo techniques to explore the viable parameter space efficiently (e.g., HYPERSPACE [26] and SloppyCell [46]). Here, we used a previously developed MCMC method [26] to explore the parameter space, initiating a random walk away from each of the 40 acceptable parameter sets and keeping 100 viable sets with an ER 10 for each point. This resulted in an ensemble of 4,000 parameter sets that could all capture our experimental observations, and allowed us to assess robustness of our predictions to parametric uncertainty. In addition, before proceeding, we ensured that all models in the ensemble satisfied our design criteria (S5 Fig).

Minimal H 2 O 2 detoxification model
Based on the tight predictions that AHP, HPI, and HPII would dominate clearance (Fig 5), we sought to determine the minimal reaction network required to capture all of our data. To do this, we adopted a previously used two-tiered approach that first deletes reactions in a random order, and then re-optimizes uncertain parameters to determine if adjusted parameters would allow deletion of additional reactions [33]. Beginning with the best model in our ensemble and using this method, we determined that 70 out of our 75 reactions could be removed without increasing the ER beyond a threshold of 10. The essential reactions to the network were the major detoxification enzymes (AHP, HPI, and HPII) and degradation of AHP and HPI. In the case of AHP, a drop in active enzyme could indicate degradation, or alternatively a decrease in available NADH, which is held constant during simulation. On the other hand, H 2 O 2 is the only substrate of catalase, which suggests that the importance of degradation reflects a decrease in concentration of functional enzyme.

Parametric sensitivity analysis
In addition to identifying reactions in the network that are dispensable to capturing the H 2 O 2 clearance data presented in Fig 5, we identified those uncertain parameters that influenced the simulations. Using the optimal parameter set, we individually varied each of the 13 optimized parameters within their bounds. Changes to the Fenton reaction rate constant and Fe 2+ and Fe 3+ initial concentrations never increased the ER to beyond 10. All other uncertain/trained parameters perturbed simulations to varying degrees, and their impact was quantified and displayed in S6 Fig

Impact of carbon deprivation on H 2 O 2 detoxification
With the model developed and an ensemble of viable parameter sets identified, we sought to assess its predictive capabilities on a physiologically-relevant environmental perturbation. There is growing evidence that microbial killing within macrophages is a combined effect of the toxic environment and a scarcity of nutrients [47][48][49]. We therefore chose to investigate how H 2 O 2 detoxification changes during carbon starvation. In the absence of an exogenous source of energy and carbon, the abundance of reducing equivalents can fall to limiting levels [50] and energetic processes such as translation can be hampered [51]. These effects could impact the stress response network by limiting AHP activity and inhibiting H 2 O 2 -dependent induction of AHP and HPI. To determine if carbon starvation in the media used here leads to NADH depletion, we directly measured NADH and NAD+ in M9 media cultures with and without glucose and found that a significant reduction in NADH occurred in carbon-starved cultures (S7 Fig). To see if carbon starvation depresses NADH to levels that inhibit enzyme activities, we measured respiration, which is an NADH-driven process, in M9 media in the presence and absence of glucose and found it to be significantly impaired when glucose was omitted (S8 Fig). In addition, to see if a lack of carbon reduces NADH to levels that impair AHP activity, we monitored clearance of 10 μM H 2 O 2 in a strain with AHP as the lone major detoxification enzyme (ΔkatE ΔkatG), and found that omission of glucose completely inhibited H 2 O 2 clearance in this strain (S9 Fig). In accordance with these results, model predictions indicated that if NADH was not held constant, AHP would drain it from the system in less than a second (S10 Fig). The impact of glucose starvation on induction of AHP and HPI expression was also assessed with the use of GFP reporter plasmids. Omission of glucose completely inhibited H 2 O 2 -dependent induction (S11 Fig). Therefore, to simulate the impact of carbon deprivation, NADH concentrations were no longer held constant and protein production was set to zero. Ensemble predictions for carbon deprivation (-glucose) were made using the complete reaction network and are shown in Fig 6A-6D, along with the carbon-replete control (+ glucose). These predictions were experimentally confirmed and the data are presented in Fig 6I-6L, orange). To quantify how the different elements of glucose deprivation (NADH limitation, inhibition of translation) contributed to the observed phenotypes, we investigated the individual effects of NADH depletion or translation inhibition with simulation controls (Fig 6E-6H). At lower treatment concentrations (10 and 25 μM H 2 O 2 ), starvation was predicted to slow detoxification as a result of NADH depletion, whereas inhibition of protein synthesis was predicted to have a negligible effect. At 100 μM H 2 O 2 , glucose-deprived cultures were predicted to clear H 2 O 2 comparably to glucose-fed cultures, with neither reducing equivalent availability nor enzyme production substantially hindering detoxification. The impact of starvation at 400 μM H 2 O 2 was predicted to be largely mediated by translation. We note that although selective inhibition of NADH production and usage was not feasible due to the wide variety of sources and sinks, targeted inhibition of translation was experimentally tractable with the use chloramphenicol (CAM) (S11 Fig). Experimental confirmation of clearance by CAM-treated cultures is shown in Fig 6I-6K. Unfortunately, CAM-treatment led to cell death at 400 μM H 2 O 2 (S1X Fig), which prevented direct confirmation of the prediction at that concentration. Interestingly, this cell death suggested that translation of some protein other than AHP or HPI is important to survival at 400 μM H 2 O 2 , because we demonstrated that carbon deprivation inhibited induction of AHP and HPI at 400 μM H 2 O 2 (S11 Fig), and it is known that carbon deprivation can have promoter-specific effects [51] and CAM stops synthesis of all proteins.

Discussion
The toxic nature of H 2 O 2 makes it an ideal weapon in inter-species warfare, and it is used as such by the immune system during infection [2] and even by other bacteria in niche competitions [15]. Bacteria have evolved numerous defense systems, which can differ significantly in their substrate requirements, reaction mechanisms, and regulation, and the complexity of these defense networks and the broad reactivity of H 2 O 2 necessitate the use of computational modeling for quantitative interpretation and prediction of H 2 O 2 distributions in cells [52]. Due to its importance as a signaling molecule, models of H 2 O 2 metabolism in mammalian systems have been constructed, and they have included enzymatic detoxification of H 2 O 2 [22][23][24], oxidation of cysteine residues [22], transport of H 2 O 2 across membranes [22][23][24], and oxidation of targets involved in signaling [24]. However, beyond their specificity for mammalian systems, none have accounted for uncertainty in optimized parameters or included synthesis or inactivation of enzymes, side reactions of H 2 O 2 , or other reactive oxygen species present in the network. In bacteria, modeling efforts have focused on subsystems affected by H 2 O 2 stress, such as that of the thioredoxin system in E. coli, which included the oxidation of thioredoxin and the reduction of oxidized thioredoxin, methionine sulfoxide, protein disulfides, and 3'-phosphoadenosine-5'-phosphosulfate [25]. These previous efforts inspired us to construct a quantitative, systems-level model of H 2 O 2 stress in E. coli that includes media and cellular compartmentspecific species and reactions; H 2 O 2 -dependent transcriptional regulation, inactivation, and activity of H 2 O 2 detoxification enzymes; reductases to reduce oxidized species; O 2 − • and •OH and their related reactions (e.g., oxidation of all twenty amino acids by •OH); and reactions of H 2 O 2 with other metabolites such as glutathione and the α-keto acid pyruvate. In addition, we addressed structural uncertainty in the model using an iterative computational and experimental methodology, and assessed parametric uncertainty using an MCMC procedure, which enabled the robustness of model predictions to be assessed. Similar ensemble approaches have become popular methods to account for parametric uncertainty [53][54][55][56][57][58][59][60][61][62], and several techniques have been developed to efficiently explore parameter spaces [26,46].
One power of quantitative computational modeling is its ability to predict emergent systems behavior [33,44,[63][64][65]. For instance, Schaber and colleagues used an ensemble of possible models describing different hypotheses regarding the mechanism of the high osmolarity glycerol (HOG) pathway in yeast to uncover novel features of the pathway [44]. Here, we leveraged our model to gain a quantitative understanding of how carbon deprivation, which bacteria encounter in phagocytes [47][48][49], affects H 2 O 2 detoxification by E. coli. Accounting for the NADH limitation and translational inhibition that occurs with carbon source starvation, our simulations were able to correctly predict H 2 O 2 clearance dynamics. Upon dissection of simulation results, delayed clearance at lower concentrations was attributed to reduced AHP activity from NADH depletion, whereas at higher concentrations carbon-starved cultures resembled carbon-replete cultures because pre-expressed HPI dominated H 2 O 2 detoxification, suggesting that both NADH availability and induction of AHP and HPI synthesis at H 2 O 2 concentrations > 25μM were of minor importance. We note that changes in concentrations of other metabolites, such as ATP, occur in carbon-starved cultures [66], and that they were not explicitly accounted for here because they did not directly act as a substrate in any of the reactions of the model. Rather we anticipate that some of the metabolite perturbations were implicitly accounted for because they contributed to the inhibition of translation and/or depletion of mean. Results for clearance of 400 μM H 2 O 2 by CAM-treated cultures were omitted based on significant cell death (S1X Fig). In all simulations, windows represent the maximum and minimum of the predictions from the 4,000 models in the ensemble. Solid lines within the window show the most likely model. NADH that were included. These data demonstrate that the nutritional status of the environment can have a major impact on bacterial H 2 O 2 defenses, but the extent of that impact depends on the quantitative level of H 2 O 2 . Beyond carbon deprivation, it would be interesting to see how other types of starvation (e.g., sulfate, iron) influence H 2 O 2 detoxification, since bacteria are subjected to oxidative stress in various scenarios [1,15,67], and we expect that dependencies distinct from those of carbon source starvation could be observed if other types of limitation influence NADH availability and translation differently.
In immune cells, phagocytized bacteria are exposed to ROS with an oxidative "burst" from NADPH oxidase, which then tapers over time [68][69][70]. In this work, we examined detoxification of a burst of H 2 O 2 with bolus treatments, and note that more complex treatment dynamics could be handled by the model developed here. For instance, the model could be adjusted for continuous treatment by adding an H 2 O 2 delivery reaction, which could be achieved experimentally with a fed-batch reactor. Alternatively, H 2 O 2 could be provided through indirect means, such as with exposure to redox-cycling agents like paraquat [71]; though obtaining accurate estimates of H 2 O 2 production from such compounds could be a challenge, because generation would be cell-dependent. Different H 2 O 2 delivery dynamics could have a profound impact on the kinetic competition for H 2 O 2 , and as long as H 2 O 2 influx can be accurately accounted for the model developed could prove invaluable for interrogating its distribution.
One area of growth for the platform we developed is adaptation of the model to allow for analysis of lethal H 2 O 2 concentrations. In its current form, the size of the cellular compartment is fixed and enzymatic reactions do not occur in the media compartment. To model lethal concentrations of H 2 O 2 , cell lysis has to be accounted for in terms of reduction in the volume and surface area of the cellular compartment and the addition of certain enzymatic activities to the media compartment, such as catalase, which can function when released from cells. In addition, it might be necessary to diversify the cellular compartment if all cells that die do not lyse but contain compromised translational and/or catabolic activities. Despite this added complexity, the ability to accurately simulate lethal damage, such as that involving DNA and the membrane, could provide insight into H 2 O 2 -induced death.
Models akin to the one described here will improve understanding of bacterial defenses against host immune responses, and possibly suggest targets for novel anti-virulence therapies [52]. For example, an existing model of nitric oxide (NO•) defenses in E. coli [33] has provided valuable predictions regarding NO• delivery rates that maximize antimicrobial activity [63]. Additionally, it provided a framework that allowed for a model-guided investigation of the underlying mechanism of an NO•-sensitizing mutation, ΔclpP, in E. coli [72]. The success of the NO• model provided inspiration to develop a similar model for H 2 O 2 , which is another toxic, diffusible metabolite used by immune cells when fighting infection [1,2]. We anticipate that the H 2 O 2 model developed here will yield novel quantitative insight into the kinetic competition for H 2 O 2 in E. coli and provide a framework for the mechanistic investigation of perturbations that affect clearance, while illuminating targets to sensitize bacteria to immune attack.

Materials and Methods
Bacterial strains E. coli K-12 MG1655 was used in all experiments. ΔkatE and ΔkatG mutations were transduced into MG1655 from their respective strains in the Keio collection [73] by the P1 phage method. The ΔahpCF mutation was generously provided by Michael Kohanski and transduced into MG1655 using the P1 phage method. All antibiotic markers were cured out using pCP20 [74]. The known antioxidant pyruvate (25 mM) was included in all LB agar plates to prevent a toxic build-up of H 2 O 2 [75]. Deletions were PCR checked for proper chromosomal integration with a forward primer external to the gene and reverse primer within the kanamycin resistance cassette (kan R ) before curing. Internal primers were used to check for gene duplication. In cured strains, external primers before and after the gene were used to check for proper scar size. All PCR primers are listed in S4 Table.

Measurement of H 2 O 2 clearance
Overnight cultures were inoculated from −80°C stocks and grown for 20 hours in 1 mL LB + 30-75 U/mL catalase (bovine liver catalase at 2,000-5,000 units/mg protein: Sigma Aldrich), then used to inoculate 20 mL M9 10 mM glucose medium + 30-75 U/mL catalase to an OD600 of 0.01 in 250 mL baffled flasks. Catalase was added to prevent the possibility of H 2 O 2 accumulation in strains lacking major detoxifying enzymes, and added to wild-type cultures to maintain consistency across strains. The catalase concentration was chosen based on the amount required to maintain growth in a mutant lacking all major detoxification enzymes (ΔkatE ΔkatG ΔahpCF), beyond which increasing the catalase concentration no longer increased growth rate or terminal cell density in the case of overnights. Cultures were grown at 37°C with shaking at 250 rpm for 8 h (OD 600 0.3-0.6, 0.15 for the slower growing ΔkatE ΔkatG ΔahpCF). After the 8 hour growth period, 12 mL of culture was removed to a prewarmed 15 mL Falcon tube and centrifuged at 37°C and 4,000 rpm for 10 min. 10.8 mL of spent media was removed, the cell pellet resuspended, and 1 mL transferred to a warm 1.5 mL microcentrifuge tube. Cells were washed a total of four times to remove all catalase. Washes consisted of spinning down at 14,000 rpm for 2 min, removing 980 μL of media, and resuspending the cell pellet with 980 μL fresh warm media. For samples lacking glucose during challenge with H 2 O 2 , glucose was omitted during the final wash step and in the inoculated flask. For CAM treatment assays, all wash steps were performed with 100 μg/mL CAM.
Prior to inoculation with washed cells, 20 mL fresh M9 10 mM glucose media in 250 mL baffled flasks were warmed to 37°C. A bolus of 10, 25, 100, or 400 μM H 2 O 2 was added to the flasks, and the time 0 point was measured, after which flasks were inoculated to an OD 600 of 0.01. At desired time points, 200 μL was removed to a 1.5 mL microcentrifuge tube and centrifuged at 15,000 for 3 min. 150 μL of the supernatant was moved to a sterile microcentrifuge tube and stored at 4°C until H 2 O 2 concentration could be measured. Samples were assayed for H 2 O 2 within 2 h of harvesting. H 2 O 2 in the supernatant was measured using the Amplex Red Hydrogen Peroxide/Peroxidase kit (Life Technologies) according to the manufacturer's instructions after dilution to below 10 μM H 2 O 2 . A standard curve spanning 0 to 10 μM H 2 O 2 was used to calculate H 2 O 2 concentrations. A fresh standard curve was produced for each Amplex Red assay to account for increasing background fluorescence over the course of the day due to the sensitivity of Amplex Red to both light and air [76].
To assess whether centrifuging to remove cells was sufficient, we compared values from this method to sterile filtering (0.22 μm pore size) of samples, as well as centrifuging + filtering, for the 30 min point in the 400 μM H 2 O 2 clearance assay (S12 Fig). Centrifuging alone was not significantly different from filtering (p = 0.45) or centrifuging + filtering (p = 0.40) based on a two-sample t-test with unequal variance. To determine if the exogenous catalase that was added to the pre-culture steps affected clearance profiles, we performed identical experiments for wild-type propagated without catalase in all steps (S13 Fig). The presence of catalase in the pre-culture did not affect wild-type clearance of H 2 O 2 at any concentration.

Measurement of colony forming units (CFUs)
To determine whether H 2 O 2 treatment resulted in cell death, we quantified CFUs throughout the clearance assays. After isolating the H 2 O 2 -containing supernatant for Amplex Red assays as described above, an additional 30 μL supernatant was removed from centrifuged samples and discarded, to achieve a greater fold-dilution of H 2 O 2 during the first wash step. In the first wash step, 980 μL of PBS was added and the cell pellet was resuspended. Samples were centrifuged again at 15,000 rpm for 3 min, 980 μL of the supernatant was removed, and the cell pellet was resuspended a final time in 80 μL PBS. Plating was performed using the serial dilution method, and samples were plated on LB agar containing 25 mM pyruvate to scavenge any residual H 2 O 2 remaining in the pellet and any endogenously produced H 2 O 2 in scavengingdeficient strains. Plates were incubated for 16 h at 37°C prior to counting colonies.

Oxygen measurement
MG1655 cultures were grown and washed identically to the H 2 O 2 clearance assay. After the final wash, the resuspended cells were used to inoculate 10 mL of pre-warmed M9 with or without glucose in a 50 mL Falcon tube containing a sterile magnetic stirring bar, immersed in a stirred water bath at 37°C, to an OD 600 of 0.1. Cells were allowed to consume oxygen for ten minutes before being treated with 5 mM KCN to halt respiration, which consumes the majority of O 2 in E. coli cell cultures under these conditions [63]. The percent oxygen saturation was measured at a frequency of one reading per second using the FireStingO2 fiber-optic O2 meter with the OXROB10-CL2 robust oxygen miniprobe (PyroScience, GmbH). Temperature fluctuations were compensated for using the TDIP15 temperature sensor (PyroScience GmbH) and the FireSting Logger Software. The equilibrium oxygen concentration was used to convert the percent saturation to concentration, and was determined by calibrating the probe in ultrapure Milli-Q water at 37°C, which has an oxygen concentration of 210 μM [77], and transferring the probe to air-saturated M9 media. The equilibrium concentration of the media matched that of the ultrapure water.

NAD+/NADH measurement
Overnight cultures were inoculated and grown identically to the H 2 O 2 clearance assay, and used to inoculate 20 mL M9 10 mM glucose medium + 30-75 U/mL catalase to an OD600 of 0.01 in 250 mL baffled flasks. Cultures were grown at 37°C with shaking at 250 rpm to an OD 600 of 0.2 (~6.5 h). Four 1 mL aliquots were transferred from the flask to warm, 1.5 mL microcentrifuge tubes and centrifuged at 15,000 rpm for 3 min. The media was removed, and the pellets were resuspended with 1 mL fresh M9 with or without 10 mM glucose and transferred to warm test tubes. The tubes were then incubated at 37°C with shaking at 250 rpm for 60 min. The time 0− point was taken directly from the flask prior to centrifuging and resuspension. NAD+ and NADH were measured using the EnzyChrom NAD/NADH Assay Kit (BioAssay Systems) following the manufacturer's protocol, except for a brief sonication step. For each measurement, 400 μL (NAD+) or 800 μL (NADH) of cell culture was transferred from the flask (time 0) or test tube (60 min) to a 1.5 mL microcentrifuge tube. The tubes were centrifuged for 3 min at 15,000 rpm, and 380 μL (NAD+) or 780 μL (NADH) of supernatant was removed and discarded. The cell pellets were resuspended in 100 μL of either NAD or NADH extraction buffer and sonicated for 20 s at room temperature at an amplitude of 10 using a Fisher Scientific Model 50 Sonic Dismembrator. The extracts were heated at 60°C for 5 min, before adding 20 μL of assay buffer and 100 μL of the opposite extraction buffer. Samples were then vortexed briefly and centrifuged for 5 min at 14,000 rpm. The supernatants were used for the NAD/NADH assay following manufacturer's instructions. An NAD+ standard curve from 0-2 μM was generated each day and used to convert absorbance to concentration. The standard curve underwent extraction protocols identical to cell samples, including sonication and heating. Since NAD+ and NADH produce identical standard curves, and NADH is more unstable, only NAD+ was provided in the kit and was used to convert both NAD+ and NADH absorbance to concentration.

Model framework
The modeling framework used in this work largely followed that used by Robinson and Brynildsen [33]. It is composed of a system of ordinary differential equations that are numerically integrated to provide predicted species concentration over time. We begin with a mole balance for all species in the model: where N is an s x 1 vector representing the total amount of given species in moles, S is the s x r stoichiometric matrix, and v is an r x 1 vector representing reaction rates in moles per time.
Here, s indicates the number of species in the model, and r indicates the number of reactions in the model. Since most reaction rates are calculated on the basis of concentration per time, the rate vector was converted into these units in the following manner.
where V rxn is an r x r diagonal matrix representing the volumes of the compartments in which the reactions are taking place: V cell for intracellular reactions, V media for reactions taking place only in the media, and V total for exchange reactions. Due to experiments also being performed on a concentration basis, N was converted into units of concentration by performing the following operation.
where V spec is a diagonal s x s matrix of species compartment volumes: V cell for intracellular species, V media for species in the media compartment, and V total for species that freely diffuse across the cell membrane (e.g., O 2 , H 2 O 2 in non-gradient models). The left-hand side of the equation is equivalent to dC/dt when the volume does not vary appreciably over the course of the experiment. To avoid having a culture-volume-specific model, we transformed the volume dependencies into volume fractions by multiplying and dividing by V total (V cell /V total , V media / V total , V total /V total ) and rearranging the equation with the use of the commutative property of scalar multiplication.
where F spec is an s x s diagonal matrix of the volume fractions for species and F rxn is an r x r diagonal matrix of the volume fractions for reactions. By making this adjustment, we avoid the need of requiring total culture volume as an input into the model, and simplify the input to optical density (OD 600 ) that can readily be converted to volume fractions.

Initial species concentrations
Most initial species concentrations were obtained from literature (S1 Table) [21,25,[79][80][81][82][83][84][85][86][87][88][89][90]. The equilibrium concentration of oxygen was determined by calibrating a FireStingO2 fiberoptic O 2 meter with the OXROB10-CL2 robust oxygen miniprobe (PyroScience, GmbH) in ultrapure Milli-Q water at 37°C and transferring it to air-saturated M9 10 mM glucose medium at 37°C. We calculated the concentration in our media by comparing it to the known value in deionized water [77], which is 210 μM. The value in our media was equivalent. The initial H 2 O 2 concentration was set to the initial average value of the data when optimizing parameters (e.g., 25.89 μM instead of 25 μM for wild-type) to avoid penalizing model fit for experimental error. When making forward predictions, the concentration was set to the anticipated initial concentration (e.g., exactly 25 μM). Initial species that were trained on experimental data included AHP, HPI, HPII, Fe 2+ , and Fe 3+ . While experimental measurements on AHP [21], HPI [82], and HPII [82] are available, their concentrations vary with environment and growth phase, as shown in our reporter assays (S11 Fig). AHP, HPI, and HPII are the major H 2 O 2 detoxification systems in E. coli. We therefore allowed flexibility in their initial concentrations, constraining them to be within the range of 0-20 μM. Because individual concentrations of Fe 2+ and Fe 3+ are unresolved, but experimentally measured to have a combined concentration of 10 μM, we allowed their initial concentrations to vary from 0 to 10 μM.

HPI and HPII reaction mechanism
Catalase activity follows a ping-pong mechanism, reacting with one H 2 O 2 to form a reactive intermediate, followed by reaction with a second H 2 O 2 molecule to return the enzyme to its original form [91,92]. Given that the substrate in the first and second reaction is the same, the rate equation simplifies to a Michaelis-Menten type structure. While inhibition of catalase activity by H 2 O 2 becomes apparent at high concentrations (e.g., greater than 100 mM H 2 O 2 for E. coli HPII [92]), it is assumed negligible at the concentrations used in our experiments, and Michaelis-Menten kinetics are appropriate. Rate equations and constants can be found in S3 Table (Reactions 69 and 70). While HPI has the ability to utilize other reducing agents at low H 2 O 2 concentrations, this activity is significantly slower than its catalase activity (about 1% the k cat of its catalase activity [21]), so the peroxidase activity was assumed to negligibly contribute to H 2 O 2 detoxification in this study.

AHP reaction mechanism
The AHP reaction cycle begins when the peroxidatic cysteine of AhpC reacts with H 2 O 2 to form a sulfenic acid, which resolves to form a disulfide bond with another cysteine residue. The active AhpC is regenerated by its reductase partner AhpF, which uses NADH as an electron donor [38]. We modeled this cycle using ping-pong reaction kinetics, with H 2 O 2 and NADH as the substrates, a structure which has been used previously [38]. Kinetic parameters were not available for E. coli, but a protein BLAST search [93] revealed 98% protein sequence identity for AhpC and 95% identity for AhpF between E. coli MG1655 and S. Typhi, so available parameters for S. Typhi were used [38,94]. Additional information and rate constants can be found in S3 Table (Reaction 71).

Expression of detoxification enzymes
In this work, the main experimental variable was the concentration of H 2 O 2 , and therefore, we opted for simplicity and only H 2 O 2 -dependent regulation of gene expression was considered. The expression of catalase HPI and AHP increase in response to H 2 O 2 [17,95], whereas the expression of catalase HPII is not dependent on H 2 O 2 [17,18]. These dependencies were confirmed using transcriptional reporters for ahpC, katE, and katG, and measuring fluorescence on a per cell basis using flow cytometry (S11 Fig). Following previous dynamic models, gene expression was modeled using a Hill equation with a coefficient of n = 1 [33,34], except for HPII which had initial concentration but was not expressed further. In addition, transcription was assumed to be limiting in the production of active enzyme, as assumed previously [33,34], and the bioavailability of ferroheme b, which is an essential cofactor of HPI, was assumed to not be rate limiting. HPI and AHP are expressed according to Reactions 67 and 68, respectively, in S3 Table. The maximum expression rate and Hill equation constants K AHP-exp,H2O2 and K HPI-exp,H2O2 are informed during parameter optimization. The bounds on the maximum expression rates are based on the highest and lowest maximum expression rates found in the work of Kotte and colleagues [34], and have been used previously when optimizing unknown expression rates [33]. Bounds on K AHP-exp,H2O2 and K HPI-exp,H2O2 were approximated by the work of Kotte and colleagues [34], which varied from approximately 2 nM to 1 mM. Here, we allowed variation from 0 to 1 mM. We note that while the parsimonious treatment of H 2 O 2dependent expression was sufficient in this work, exploration of new environments could necessitate a more comprehensive modeling of transcriptional regulation, since expression of the detoxification enzymes in this work are known to depend on numerous regulators (e.g., OxyR, RpoS, FIS, Fur) [96].

Modeling of enzyme degradation
In previous studies, enzymes have typically been modeled as undergoing first order degradation with a universal constant [33,34]. However, rates can vary greatly among different proteins [35,36], and some evidence suggests that H 2 O 2 detoxification enzymes, such as catalase and alkyl hydroperoxidase, can be poisoned by their own substrate [37][38][39]. Inactivation of catalase has been described using bimolecular kinetics [40], but it has also been suggested that poisoning should be of the same general form as the enzyme's reaction kinetics [37]. The AhpC component of alkyl hydroperoxidase is reduced by AhpF after oxidation by H 2 O 2 . If there is not enough NADH present or AhpC reacts with more H 2 O 2 before encountering AhpF, the cysteine sulfenic acid formed by the first oxidation can be further oxidized to sulfinic or sulfonic acid, rendering it inactive [38,97]. Whether AhpC inactivation is significant at the concentrations used in this work was uncertain [39]. We therefore allowed it to be degraded in a first order or bimolecular manner. Due to the indeterminacy of how these enzymes are degraded, we included all of these possible degradation scenarios. Rate equations can be found in S2 Table. We note that since AHP requires a co-substrate, NADH, and that we assume that NADH is constant (unless otherwise noted), if AHP inactivation were found to be important it could reflect degradation of AHP or reduced availability of NADH (Reaction 71).
Parameters or bounds on parameters for enzyme degradation/deactivation rate equations varied with the method of degradation. For constant degradation with a constrained constant, we used the "general" protein degradation rate reported by Kotte and colleagues [34]. When optimized, the constant degradation rate had bounds set by the longest [35] and shortest [36] protein half-life we found in literature. For both bimolecular and the more complex inactivation, bounds were based loosely on rate constants found for Aspergillus niger and bovine catalase [37,40]. Based on the gross difference between the organisms tested in those studies (fungus and mammal) and our own (bacteria), as well as the orders of magnitude difference between the Aspergillus niger and bovine catalase rates in both the bimolecular and more complex kinetics studies, these parameters were constrained to two orders of magnitude lower than the lowest reported value, and two orders of magnitude higher than the highest reported value.

Modeling of H 2 O 2 gradient across the membrane
While H 2 O 2 rapidly diffuses across bacterial membranes at a rate similar to that of water, there is evidence for and against the existence of a gradient across the membrane [41,42]. For example, an Ahp−Kat+ strain cocultured with Ahp−Kat− in the presence of a low H 2 O 2 concentration can outcompete its scavenging deficient neighbors and multiply under the stress, whereas the deficient strain only grows after the catalase proficient strain has cleared the environmental H 2 O 2 [41]. On the other hand, dilute suspensions of catalase proficient strains are readily killed by high concentrations of H 2 O 2 similarly to scavenging deficient strains, while high-density catalase proficient strains can not only survive challenge with H 2 O 2 but also protect deficient neighbors [42]. There have been strides in the ability to measure H 2 O 2 intracellularly, with the introduction of genetically-encoded indicators (HyPer) [98] and the ability to use Amplex Red intracellularly with expression of a mutated ascorbate peroxidase [99]. However, the dependence of HyPer fluorescence on reductase activity, the impact of HyPer on cellular scavenging capacity [100], and the difficulties associated with converting measurements from either method to absolute H 2 O 2 concentrations led us to use measurements of external H 2 O 2 and statistical metrics (AIC) to assess the suitability of modeling the system with or without a gradient. In one set of models, the intracellular and extracellular H 2 O 2 concentrations were equal. In the other set, we allowed for a gradient by modeling transport across the membrane as a convective mass transport process, with the effective mass transport coefficient being an additional parameter for optimization. The lower bound on the effective mass transfer coefficient was set as the permeability coefficient of H 2 O 2 across E. coli cell membranes in unstirred culture [41], adjusted for cell area and the cell density for our system. The upper bound was set two orders of magnitude higher than the permeability, to account for increased mass transfer in our chaotic shake flask system.

Parameter optimization
Parameters were optimized by minimizing the SSR using the built-in MATLAB function lsqcurvefit. Because experimental error varied for each time point, we weighted each data point's contribution to the SSR by the inverse of the variance of that point [43,[101][102][103]. The initial concentration of H 2 O 2 in the model was changed to match the experimental data before optimizing. Due to the nonlinear nature of the optimization, each model structure was initialized with random parameter sets within the defined bounds a total of 1,000 times. The progression of minimum SSR found by each of these iterations is shown in S14 Fig. A plateau suggests that additional iterations would not lead to substantial improvement in the fit. The slowest converging optimization reached a plateau after 639 optimizations.
The number of parameters optimized varied according to model structure, as presented in Table 1. Parameters related to HPI and AHP expression (4 parameters total) were optimized in all structures. The Fenton reaction rate constant varied in literature (1 parameter), and intracellular Fe 2+ and Fe 3+ concentrations (2 parameters) were unresolved. Additionally, initial concentrations of the major detoxifying enzymes (3 parameters) were unknown and can vary with growth environment and stage of cell growth. Beyond these 10 parameters, all structures that did not have constant enzyme degradation with a universal degradation constant added 3 parameters, and including an H 2 O 2 gradient added 1 parameter (a convective mass transfer coefficient).

Model ranking and selection
The introduction of additional parameters, such as enzyme degradation rate constants and mass transport coefficients, has the potential to improve fit solely by increasing the flexibility of the model. To account for the utility of additional parameters, we ranked models based on their evidence ratios (ER), or likeliness relative to the most likely model in the set. For each model, we calculated its Akaike Information Criterion corrected for small sample size (AIC c ) [43]: where n represents the sample size and K is the number of estimable parameters. Here, n is the number of data points used in the fitting procedure, and K is the number of model parameters plus 1 because regression estimates SSR and parameter values [43]. We account for unequal variances within the data by using weighted least squares, where each point is weighted by the inverse of its variance. The weight of evidence for a given model in a set of M models is given by the following [43]: where Δ i = AIC i -min(AIC). With this, an ER can be calculated, which represents the relative likelihood of a model compared to the best model in the set [27]: A larger ER indicates a more unlikely model. In this work, models with an ER greater than 10 were discarded, a cutoff that has been used previously to discard models during model selection [29].

Generation of the ensemble
To account for parametric uncertainty when making predictions, we generated an ensemble of parameter sets that all predicted the data within an ER 10. We initially attempted to use the software HYPERSPACE [26], which is a three-step process that provides a uniform sampling of the viable parameter space. However, the calculation of our cost function was computationally expensive, which lengthened the time required per iteration, and the method had not converged within 100,000 iterations. Therefore, we utilized the pre-existing Markov chain Monte Carlo function within the HYPERSPACE software, but started it from all 40 viable parameter sets that met our design criteria. We allowed approximately 200 random steps away from each viable point, and randomly selected 100 of those parameters sets with an ER 10 from each random walk. This process generated 4,000 parameter sets that had an ER 10.

Identification of the minimal H 2 O 2 model
We identified the minimal model that was able to capture our data (ER 10) using a previously developed two tiered approach [33]. In the first tier, reactions were removed from the best model in a random order and the ER was calculated. If the deletion of a reaction increased the ER above its threshold of 10, it was returned to the model and the process continued through the remaining reactions. This random deletion process was repeated 100 times. In the second tier, parameters were re-optimized after the deletion of each of the remaining reactions. If the optimization produced a model that returned below an ER of 10, the parameters were changed and the process continued. Reaction fluxes for the 35 acceptable models after fitting on wild-type data (Fig 3). E-H. Reaction fluxes for the 965 acceptable models after fitting simultaneously on wild-type and ΔkatE ΔkatG data (Fig 4). I-L. Reaction fluxes for the 40 acceptable models after fitting on wild-type, ΔkatE ΔkatG, ΔkatE, and ΔkatG data (Fig 5). Each line represents the prediction from a single model. Reaction fluxes for the 965 acceptable models after fitting simultaneously on wild-type and ΔkatE ΔkatG data (Fig 4). E-H. Reaction fluxes for the 40 acceptable models after fitting on wild-type, ΔkatE ΔkatG, ΔkatE, and ΔkatG data (Fig 5). Each line represents the prediction from a single model.  (Fig 4). Wide distributions on clearance dynamics suggest that these single mutants could be used to discriminate between models. (TIF) Exponentially growing cells were washed, resuspended in media with or without glucose, and used to inoculate M9 media +/glucose to an OD 600 of 0.1. A higher density than that used for H 2 O 2 clearance assays was necessary to observe a measureable drop in O 2 . Cells were allowed to consume oxygen for ten minutes before being treated with 5 mM KCN to inhibit respiration. The solid lines show the average of three biological replicates, and windows represent the standard error of the mean. We found that cultures in M9 media with 10 mM glucose efficiently consumed oxygen via respiration, whereas glucose deprived cultures consumed very little. To explore whether omitting glucose from the media effectively eliminated AHP activity in a regime where it dominates (10 μM H 2 O 2 ), we compared H 2 O 2 clearance in ΔkatE ΔkatG in M9 minimal media with glucose (blue), without glucose (orange), and with glucose and 100 μg/mL CAM (black) to ΔkatE ΔkatG ΔahpCF (red) and cell-free controls (green). Removal of the two other major detoxification systems leaves AHP, which requires one NADH for every reaction cycle, as the only major detoxification system. When glucose is omitted from the media, H 2 O 2 clearance is eliminated, as evidenced by ΔkatE ΔkatG in M9 minimal media without glucose never differing significantly from ΔahpCF ΔkatE ΔkatG or a cell-free control based on a twotailed t-test with unequal variance. These results were not due solely to inhibition of translation, another effect of glucose starvation (S11 predictions of the 4,000 models in the ensemble. The solid line indicates the prediction made by the most likely model. Note that the x-axis here is in seconds, not hours. We note that in all other simulations, the NADH concentration was held constant (S1 Table). (TIF) S11 Fig. Synthesis of GFP under control of katE, katG, and ahpC promoters. Wild-type cells transformed with pUA66 P katE -gfp (A-D), pUA66 P katG -gfp (F-I), or pUA66 P ahpC -gfp (J-M) were exposed to H 2 O 2 in M9 10 mM glucose (blue), M9 lacking glucose (red), or M9 10 mM glucose + 100 μg/mL CAM (green) media. We confirmed that the lack of expression of GFP from pUA66 P katE -gfp after exposure to H 2 O 2 was not due to a defect in the vector by including a 16 h overnight control in M9 10 mM glucose. As expected, katE expression increases in stationary phase as determined by an increase in fluorescence after 16 h (0 μM H 2 O 2 panel). In all cases, an empty vector control (black) was included to account for auto-fluorescence. Two biological replicates were analyzed on different days for each experiment. Each parameter set optimization was randomly initialized 1,000 times. Presented here is the minimum SSR that had been found at each of the 1,000 iterations (model index) for the following optimizations performed: model structure 1 (A), 2 (B), 3 (C), 4 (D), 5 (E), 6 (F), 7 (G), 8 (H), 9 (I), 10 (J) trained on wild-type data; model structure 2 (K) and 3 (L) trained on wild-type and ΔkatE ΔkatG data; model structure 3 (M) on wild-type, ΔkatE ΔkatG, ΔkatE, and ΔkatG data. This data suggested that more than 1,000 initializations would provide very little return for additional computational time invested. (TIF) S1