{'fraction_non_alphanumeric': 0.03362774864292943, 'fraction_numerical': 0.015571809734106174, 'mean_word_length': 5.046870653685675, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of benzoate metabolism in Pseudomonas putida KT2440

Soil microorganisms mineralize lignin-derived aromatic carbon sources using oxidative catabolic pathways, such as the β-ketoadipate pathway. Although this aromatic pathway is one of the best-studied pathways in biochemistry, the complete pathway, including its regulation by aromatic carbon sources, has not been integrated into the metabolic network. In particular, information about the in vivo operation (e.g., kinetics and flux capacity) of the pathway is lacking. In this contribution, we use kinetic modeling and thermodynamic analysis to evaluate the in vivo operation of this key aromatic multi-step pathway. The resulting ab initio deterministic model of benzoate degradation via the β-ketoadipate (ortho-cleavage) pathway in Pseudomonas putida KT2440 is presented. The kinetic model includes mechanistic rate expressions for the enzymes and transport processes. The design and experimental validation of the model are driven by data generated from short-term perturbation experiments in a benzoate-limited continuous culture. The results of rigorous modeling of the in vivo dynamics provide strong support for flux regulation by the benzoate transporter and the enzymes forming and cleaving catechol. Revisiting the β-ketoadipate pathway might be valuable for applications in different fields, such as biochemistry and metabolic engineering, that use lignin monomers as a carbon source.


Introduction
In nature, the metabolic conversion of lignin to aromatic monomers occurs via a series of cometabolic and secondary metabolic events that involve both fungi and bacteria (Bugg et al., 2011). Understanding the innate metabolism of microbes while they utilize complex aromatic compounds as their sole source of carbon and energy may provide clues for decoding the possible links among lignin degradation, aromatic monomers, central carbon metabolism and bio-product formation. Soil microorganisms (e.g., the Pseudomonads) are known to exhibit a common procession that is termed the 'metabolic funnel' to mineralize structurally diverse aromatic molecules into a few central aromatic intermediates through peripheral catabolic pathways. Those aromatic intermediates are further processed via the central aromatic pathways to the core metabolism of the cell (Díaz, 2004). The genomic analysis of a well-characterized member of the metabolically versatile bacterial species Pseudomonas putida revealed the existence of at least six major aromatic pathways for degrading central aromatic intermediates (i.e., protocatechuate, catechol, homogentisate, phenylacetate, nicotinic acid and gallic acid) (Jiménez et al., 2008(Jiménez et al., , 2002Nogales et al., 2011). Of these major pathways, the pathway involved in catechol degradation (i.e., the β-ketoadipate pathway) is the most prevalent and extensively studied. In addition to having a wide taxonomic distribution among both the prokaryotic and eukaryotic groups, the β-ketoadipate pathway provides a mechanism for the dissimilation of aromatic and hydroaromatic compounds (Parke et al., 2000). The pathway follows a two-step strategy: a ring cleavage by ring modification reactions is followed by ring fission and subsequent reactions that lead to the generation of TCA cycle intermediates (Harwood and Parales, 1996). This versatile metabolism of aromatic compounds has now sparked interest in renewable resource-based chemistry, and the possible impacts were reviewed previously (Fuchs et al., 2011;Wells and Ragauskas, 2012). One example is the microbial synthesis of cis,cis-muconic acid from benzoate. This route is explored as an alternative to chemical synthesis for the production of adipic acid, which is a commodity chemical for the production of nylon-6,6 (Mizuno et al., 1988; van Duuren et al., 2011;Vardon et al., 2015). The potential to apply these catabolic pathways for aromatic compounds in lignin- enriched streams and lignin-derived aromatic monomers can be further explored by examining the formation of medium-chain (C 6 -C 14 ) polyhydroxyalkanoates by the versatile soil bacterium P. putida KT2440 (Linger et al., 2014).
In P. putida KT2440, the β-ketoadipate pathway is native, encoded on the chromosome, and lacks the meta-cleavage route that is localized on a mega plasmid (Nelson et al., 2002). Previous studies have shown direct evidence for the activation of the orthocleavage route of the β-ketoadipate pathway during the growth of P. putida on benzoate as the sole source of carbon and energy (Cao and Loh, 2008;Ornston and Stanier, 1966). Although numerous biochemical studies of the in vitro properties of the corresponding enzymes are available, the in vivo kinetics of this pathway remains poorly understood.
The central goal of this study is to quantify the in vivo metabolic response during benzoate degradation via the β-ketoadipate (ortho-cleavage) pathway in P. putida KT2440 (Sudarsan et al., 2014). Consequently, we analyze the enzyme kinetics in vivo using a mathematical model based on deterministic rate equations. The analysis is data-driven, and the response of the metabolome to dynamic system excitation  is used to identify the kinetics of multistep enzyme reactions and to estimate the kinetic parameters using a bottom-up approach (Bruggeman and Westerhoff, 2007). The identified model is then used to predict system dynamics (i.e., changes in the intra-and extracellular concentrations of metabolites both at steady state and after a benzoate pulse). In the second part of this study, we draw on quantitative metabolome analysis and apply network-embedded thermodynamic (NET) analysis (Zamboni et al., 2008) as a tool to provide novel insights about the active regulatory sites of the βketoadipate (ortho-cleavage) pathway coupled to central carbon metabolism. The systematic incorporation of insights from the identified dynamic model and thermodynamic analysis will be useful in the future as a framework for formulating large-scale kinetic models.

Material and methods
2.1. Strain, medium, and pre-culture conditions P. putida strain KT2440, which is a plasmid-free derivative of P. putida mt-2 (Bagdasarian et al., 1981), was used throughout this study. Cells from glycerol stocks were recovered on LB agar plates. Liquid cultures were started on LB medium, followed by two subsequent cultivations on modified M9 minimal medium supplemented with benzoate. Erlenmeyer flasks (1 L) were used in a rotary shaker at 30°C and 150 rpm, ensuring fully aerated growth conditions.
The composition of the modified M9 medium supplemented with benzoate for chemostat cultivations was determined based on elemental balancing for the production of 5 g dry weight biomass per liter, with the following composition per liter: 15 g of Na-Benzoate; 2 g of Na 2 SO 4 Á 10h 2 O; 2.68 g of (NH 4 ) 2 SO 4 ; 1 g of NH 4 Cl;

Start-up of continuous cultivation and criteria for steady state
Cultivation was initiated in batch mode at 30°C in a stirred tank bioreactor with a working volume of 1.5 L the airflow was maintained at 2 L min À 1 , and the pH was maintained at 7.0 by the addition of 20% (v/v) H 3 PO 4 . Batch cultivation was monitored continuously through changes in the dissolved oxygen (DO) concentration and increases in the OD 450 . Fed-batch operation was initiated with exponential feeding of benzoate solution (100 g L À 1 ) after the complete consumption of benzoate in the growth medium (marked by a steep increase in the DO concentration). Exponential feeding was performed at a growth rate of 0.65 h À 1 and was controlled by computer-aided feed-forward control. After the desired biomass concentration of approximately 5 g CDW L À 1 (corresponding to an OD 450 of 30) was reached, the continuous mode of operation was started with a dilution rate of 0.1 h À 1 , which was increased to 0.2 h À 1 . Steady-state conditions were confirmed when the biomass concentration and the exhaust analysis remained constant over a period of five mean residence times.
For the determination of the biomass concentration, 10 mL of the cell suspension was centrifuged (Heraeus Megafuge 1.0R, rotor 8030, 15 min, 3220g, 4°C) in pre-weighed glass tubes. The biomass pellet was washed twice with a 0.9% NaCl solution, dried at 105°C until reaching a constant weight and then weighed.

Benzoate pulse experiments
To avoid the adaptation of P. putida KT2440 to benzoate-limited conditions, stimulus-response experiments were performed for a low and defined number of generations. The steady-state conditions at a dilution rate of 0.2 70.003 h À 1 were perturbed by increasing the residual concentrations of benzoate with a defined pulse of benzoate (introduced manually with a syringe) to generate a bulk concentration of 1.3 g L À 1 . Samples were withdrawn rapidly from the bioreactor at regular intervals. The opening times of the valves were synchronized and controlled via a computeraided automated sampling system. During sampling, the feeding rate was corrected with a computer-aided peristaltic pump to compensate for the loss of reactor volume and to maintain a constant dilution rate of 0.2 h À 1 throughout the sampling period.

Sampling procedure
Intracellular metabolites were quantified using a differential method of calculation (Taymaz-Nikerel et al., 2009;Theobald et al., 1997), in which the intracellular metabolite concentrations were calculated from the difference between the metabolite concentrations obtained in the total broth extracts and cell-free extracellular samples. The specific cell volume was determined to be 3.8 mL g CDW À 1 for P. putida KT2440 under the above-mentioned culture conditions. For the measurement of the total metabolite concentrations, 4 mL of sample broth was withdrawn, quenched rapidly in 1 mL of HClO 4 ( À 25°C, 35% w/v) and extracted for 15 min by end-overend mixing at 4°C. The compounds benzoate cis-diol and cis,cismuconate were found to be acid labile, whereas benzoate and catechol were stable under acidic conditions. Consequently, the pH of the sample was adjusted to either pH 7 (for analyzing benzoate, benzoate cis-diol and cis,cis-muconate) or an acidic pH of 2 (for analyzing catechol). After centrifugation (Heraeus Megafuge 1.0R, rotor 8030, 15 min, 3220g, 4°C), the supernatant was collected and stored at À 20°C until further analysis using HPLC.
Four milliliters of sample broth was sampled and filtered manually at room temperature to determine extracellular concentrations. The cells were separated rapidly using glass microfiber filters (Whatman, GD/X 25 syringe filter, 0.45 mm), and the obtained filtrate was stored at À 20°C until further analysis. The obtained filtrate was immediately mixed with 1 mL of HClO 4 to prevent the auto-oxidation of catechol in the culture filtrates. The exact volumes of sample processed were quantified gravimetrically by weighing the tubes at every individual step of sample processing.

Metabolite analysis
Metabolites, such as benzoate, benzoate cis-diol, cis,cis-muconate, catechol and β-ketoadipate, were analyzed and quantified using reverse phase HPLC. Except for β-ketoadipate, all other metabolites were separated on a CC Nucleosil C-18 HD column (100 Å pore size, 5 mm particle size, 25 cm Â 4 mm inner diameter; Macharey-Nagel, Oensingen, Switzerland) and detected via an HPLC diode array detector (Elite LaChrom, Merck-Hitachi, Germany), with wavelengths ranging from 200 to 400 nm. The conditions used to separate the above-mentioned metabolites are described below. Solvent A was 100 mM ortho-phosphoric acid (pH 2), and solvent B was methanol. The starting concentration was 10% B in A. After 5 min, a gradient was established by increasing the concentration of B by 3% per min at a flow rate of 1 mL min À 1 . A sample volume of 10 μL was injected, and the column temperature was kept constant at 25°C. The maximum wavelengths were 230 nm for benzoate, 262 nm for benzoate cisdiol, cis,cis-muconate and its isomers, and 276 nm for catechol. Under acidic conditions, benzoate cis-diol was found to decompose into salicylate and phenol; however, this phenomenon was not observed at neutral pH. Irrespective of handling conditions, cis, cis-muconate was found to mostly degrade spontaneously into its cis-trans isomer. Hence, the degradation kinetics of cis,cis-muconate was assessed, and the exact amount of cis,cis-muconate was quantified by combining the peak areas of cis,cis-and cis-transmuconate, with a conversion factor of 0.855 for the cis-trans isomer.
The compound β-ketoadipate was separated on a Luna C8 (2) column (100 Å pore size, 5 mm particle size, 150 mm Â 4 mm inner diameter; Phenomenex, Torrance, USA) and detected via an HPLC Corona charged aerosol detector (Dionex Softtronic GmbH, Germering, Germany). Separation was achieved using mobile phase A, which consisted of 1% acetonitrile and 0.4% tri-fluoroacetic acid in distilled H 2 O, and acetonitrile as mobile phase B. The following elution profile was employed: elution started with 2% B for the first 8 min, followed by a sharp gradient of 98%, and elution was maintained until 14 min at a flow rate of 0.8 mL min À 1 . To equilibrate the column back to its initial conditions, the gradient was decreased to 2% B after 14 min. To quantify β-ketoadipate, 20 μL of the sample was injected, and the column temperature was kept constant at 40°C.

Stoichiometric model of benzoate metabolism
A stoichiometric metabolic network that included the reactions of the β-ketoadipate (ortho-cleavage) pathway, the TCA cycle, gluconeogenesis and the pentose phosphate pathway was compiled from the biochemical knowledge generated by the genomescale metabolic model of P. putida (Nogales et al., 2008). These pathways can be regarded as four different modules that were aggregated to describe the steady-state behavior of benzoate metabolism. Although P. putida differs from Escherichia coli in terms of hexose utilization, the percent carbon content in the dry cells of both strains is approximately 47%, and the values do not differ significantly (Bratbak and Dundas, 1984). Van Duuren et al. (2013) also showed that the differences in the biomass composition of the two strains are negligible. Hence, the precursor requirements for biomass formation in P. putida were assumed to be similar to those of E. coli (Emmerling et al., 2002) in this study. In this model, we used a P/O ratio of 1, according to the work of Blank et al. (2008). In the simulation run with an alternative P/O ratio of 2, only the ATP production rate changed, while the fluxes in the benzoate pathway remained unchanged. To estimate the net fluxes, the NETTO module of the FiatFlux software (Zamboni et al., 2005) was used, and the absolute in vivo (net) fluxes were estimated using global material balances that were derived using (a) the stoichiometric reaction matrix, (b) physiological data (substrate uptake rates, product formation rates, and growth rates), and (c) the precursor requirements for biomass synthesis. Fischer et al. has described the appropriate statistical data treatment for flux estimation (Fischer et al., 2004), in which the residuals were weighted by division with the experimental variances s 2 from extracellular measurements. For each metabolite that was used as a precursor for biomass synthesis, a proportional error of 5% was assigned. The fluxes were then estimated by minimizing the sum of the weighted square residuals from both the metabolite balances and flux ratios using the MATLAB function fmincon starting from a randomly generated flux distribution.
The benzoate metabolic network is well determined without the use of flux ratios. The steady-state substrate uptake and product formation rates were calculated from the mass balance equations of the respective metabolites. The specific benzoate uptake rate at steady state (q BZT ) was then calculated using the following equation: C BEN feed is the concentration of benzoate in the feed medium; C x is the biomass concentration, and D is the dilution rate during continuous cultivation.

Kinetic model structure
To develop a dynamic model of the β-ketoadipate (ortho-cleavage) pathway, mass balance equations for the extracellular and intracellular metabolites were formulated. The derived equations were of the following form: ij j i C i denotes the concentration of metabolite i, and μ is the specific growth rate or dilution rate during continuous cultivation. v ij is the stoichiometric coefficient for the i th metabolite in reaction j, and r j is the rate of the reaction j. The structure of the metabolic model is outlined in Fig. 1, and the complete set of mass balance equations is listed in Table 1. The concentrations of the metabolites were expressed in mM, and the rates were expressed in mM min À 1 .

Kinetic rate expressions
The sequence of enzyme-catalyzed reactions involved in the breakdown and cleavage of benzoate is described by mechanistic rate expressions ( Table 2). The rate expressions were chosen based on published in vitro investigations of the respective enzymes. For cases for which published information was scarce, several kinetic expressions were tested with appropriate assumptions, and those that represented the in vivo experimental data reasonably well were chosen. The selection of appropriate kinetic equations is described first for transport systems and then for benzoate catabolism.
observed non-specificity towards benzoate. The secretion processes of intracellular metabolites, such as benzoate cis-diol, catechol, cis,cis-muconate and β-ketoadipate, were modeled by considering the transport processes that are mediated by passive diffusion or active transport. For benzoate cisdiol and catechol transport, the structural property of the model is based on the assumption of reversible influx-efflux kinetics of catechol, with a Hill coefficient to describe the extracellular concentrations of benzoate cis-diol and catechol [Eqs. (22) and (24)]. The use of a Hill coefficient of 1 (i.e., c¼1 in Eq. (24)) failed to successfully describe the experimental observations (i.e., the transport of catechol). However, introducing an empirical Hill coefficient (i.e., c4 1) made it possible to mimic the active transport of toxic catechol through the cell membrane to the external medium; in addition, an increase in the Hill coefficient led to an apparent reduction in the effective influx of catechol. Catechols are also known to damage the integrity of the membrane due to their redox cycling activity (Schweigert et al., 2001). Although catechol is less lipophilic (with a log P o/w of 1.5), it is water soluble and can partition well into the cell membrane, causing a change in membrane fluidity (de Bont, 1998). Functional analysis of the genes mucK (Williams and Shaw, 1997) and pcaT (Parke et al., 2000) revealed that their gene products are involved in the transport of cis,cis-muconate and β-ketoadipate. The diffusion of charged compounds across the cell is described by the constant field equation (Goldman, 1943;Keifer and Roos, 1981). The ionic charges of the compounds cis,cis-muconate and β-ketoadipate are À 2. For the membrane potential, a constant value of À50 mV was assumed. With the assumption that the translocation of charged compounds through a carrier was symmetrical (according to previous studies (Lancaster and Hinkle, 1977;Michels et al., 1979)), the equation was modified to include the effect of the complex formed by the diffusion ion and the carrier. The concentration of the complex (ion: carrier protein) has been reported to be similar to that for the enzyme-substrate complex of Michaelis-Menten kinetics. Under the assumption of constant pH, the effect of protons was found to be negligible; thus, we modified the proposed equations to consider the effects of cis, , 2 , where P, T and R are the fermenter's pressure (bar), the cultivation temperature (K) and universal gas constant (bar L mol À 1 K À 1 ), respectively; is the partition coefficient of O 2 between gas and liquid (Bloemen et al., 2003); and 1 is the gas-liquid transfer coefficient for O 2. Benzoate transport: Benzoate 1,2-dioxygenase: Benzoate cis-diol dehydrogenase: Benzoate cis-diol transport: Catechol transport: cis-muconate and β-ketoadipate alone [Eqs. (28) and (29)]. Depending on the extracellular pH and the transport conditions, the membrane potential might vary; hence, our assumption of a constant membrane potential of À50 mV must be verified in future experiments. The rate expressions used for the transport of cis,cis-muconate and β-ketoadipate can be extended to include the effects of proton concentrations (Andersen et al., 2009) and the membrane potential change (Page, 1987); however, this topic is outside the scope of investigation of this study.

Benzoate catabolism
The foremost reaction catalyzed by benzoate 1,2-dioxygenase (encoded by benABC) involves ring hydroxylation and the incorporation of dioxygen into the product. This reaction was modeled as an irreversible double substrate kinetic (Ge and Eltis, 2003) [Eq. (20) a For cis,cis-muconate and β-ketoadipate transport, = γ= = ∆Ψ z -2; -1.915 (encoded by benD)), we tested several rate expressions that describe the rapid increase and decrease in the intracellular benzoate cis-diol concentrations. The in vivo behavior of benzoate cis-diol could be satisfactorily explained by assuming Ping-Pong kinetics with dead-end inhibition by catechol [Eq. (21)] (Segel, 1975;Simpson et al., 1987). The formation and consumption of the cofactor NAD þ by the above two reactions is assumed to be stoichiometrically balanced; hence, the influences of these factors on the enzyme kinetics were not investigated in this study. The catA gene is known to encode the enzyme catechol 1,2-dioxygenase, which is a key enzyme responsible for the intra-diol cleavage of the aromatic ring and the conversion of catechol to cis,cis-muconate (Nakai et al., 1988). Jiménez and colleagues previously reported the presence of a second copy of the gene catA (named catA2) within the ben cluster for benzoate degradation in P. putida KT2440 (Jiménez et al., 2002) and also recently elucidated the biological role of the gene product of catA2 (CatA2) as a metabolic safety valve that reduces the intracellular levels of toxic catechol (Jimenez et al., 2014). In this study, the roles of the gene products of catA and catA2 (i.e., in catalyzing the conversion of catechol to cis,cis-muconate) were represented by a combined kinetic rate expression [Eq. (23)], which was described by a double substrate equation with allosteric activation, where oxygen acted as both the substrate and an effector (Broderick, 1999;Bugg and Winfield, 1998). Because oxygen was a second substrate for both dioxygenases, the oxygen balance was included as a state variable [Eq. (16)], and the oxygen demands of both dioxygenases were considered to describe the rate of oxygen transfer into the cells.

Estimation of maximal rates
At steady state, the deterministic form of the rate equations can be envisaged as follows: Here, estimates of the maximal rates of the enzymatic reactions (r j max ) were obtained from the steady-state fluxes (r j ) and steadystate concentrations (C j ), and an appropriate estimate of the corresponding parameter vector (P j ) was obtained from the in vitro measurements. Subsequently, with an assumption that the enzyme concentrations remained frozen at their steady-state values after a substrate pulse , the unknown maximal rates were then given by the equation:

Model implementation and analysis
The Matlab-based software packages SBtoolbox2 and SBPD toolbox were applied to integrate the system of non-linear differential equations and to estimate the model parameters (Schmidt and Jirstrand, 2006). Parameter estimation was performed using an enhanced Scatter Search algorithm integrated into the SBPD toolbox (Rodriguez-Fernandez et al., 2006). The optimization routine integrates and minimizes the difference between the measured concentrations ( C i meas ) and the model-predicted concentrations (C i pred ), represented as the minimization of a weighted distance measure, . The set of model parameters identified from the above-mentioned procedure was subjected to parameter fit analysis, in which repeated estimations were performed from randomly perturbed initial conditions. The parameters were perturbed randomly between the low and high bounds around the identified values. The results of parameter fit analysis were then used to gain statistical insights into the estimated parameters.
LINDO API library (LINDO Systems Inc., Illinois, USA). The genomescale metabolic model of P. putida (Nogales et al., 2008) was used to predict the minimum set of flux directions for benzoate. The comprehensive summary of the thermodynamic information that was reported by Zamboni et al. (2008) was extended with additional information regarding the metabolites involved in benzoate degradation, which was obtained from group contribution methods (Finley et al., 2009). Beyond the actual measurements from benzoate-limited cultivations, the concentrations of non-measured metabolites were allowed to vary between 0.001 and 10 mM in the NET analysis. The optimization procedures for determining the feasible range of metabolic network operation using the Gibbs energies of reactions and the flux balance algorithm applied to constrain the metabolic network were summarized by Kümmel et al. (2006a) (Kümmel et al., 2006a;Kümmel et al., 2006b).

Metabolic response of P. putida to benzoate
Benzoate uptake in P. putida is driven by an active transport system (Thayer and Wheelis, 1982), which has a specific uptake rate of 3.4 mmol benzoate g CDW À 1 h À 1 (or 0.5 g benzoate g CDW À 1 h À 1 ).
Benzoate degradation via the β-ketoadipate (ortho-cleavage) pathway leads to equimolar quantities of CO 2 , succinate, and acetyl-CoA (Fig. 2) (Ampe and Lindley, 1996). Succinate and acetyl-CoA serve as fuels for anabolic reactions. The carbon skeletons that form glucose-6-phosphate originate from gluconeogenesis. The enzyme oxaloacetate decarboxylase (encoded by oadA) was considered as the first enzyme unique to gluconeogenesis. Pyruvate enters a series of predominantly gluconeogenic reactions that lead to the formation of fructose-6-phosphate and glucose-6-phosphate by fructose 1,6-bisphosphatase and phosphoglucoisomerase, respectively. A prominent flux via malate dehydrogenase was evident from the high flux through the glyoxylate cycle, which could be attributed to the post-translational modification of isocitrate dehydrogenase on acetogenic carbon sources (Walsh and Koshland, 1985). Dynamic changes in the extracellular and intracellular metabolites of benzoate degradation were quantified after the residual concentrations of benzoate was increased using a defined pulse of benzoate. The extracellular metabolite fraction (i.e., the ratio of the extracellular concentration (EX) of a metabolite in the filtrate to the total concentration (T) of the same metabolite measured in PCA extracts) of the benzoate pathway metabolites was high. From this analysis, it is not possible to derive any reasonable values for intracellular concentrations due to the exponential increase in the standard deviations (Vielhauer et al., 2011). In our study, this problem applied to the compounds benzoate, cis,cis-muconate and β- ketoadipate. However, in dynamic conditions, calculations of the intracellular concentrations of benzoate cis-diol, catechol, and cis,cismuconate were feasible. A step change in benzoate availability resulted in dynamic changes in the intracellular benzoate cis-diol and catechol concentrations, while the cis,cis-muconate concentrations increased steadily. The changes that propagated in the intracellular metabolites were on the order of seconds to minutes. For the steadystate conditions, the secretion of metabolites was as follows: β-ketoadipate4benzoate cis-diol4cis,cis-muconate4catechol; the extracellular β-ketoadipate concentrations approached a maximum concentration of 2 mM. These observations stimulate discussions about a metabolic overflow caused by the limited capacity of some of the enzymatic steps of the β-ketoadipate (ortho-cleavage) pathway in P. putida KT2440.

Evaluation of the kinetic model in predicting the benzoate metabolism response
To elucidate the complex interactions of the individual metabolites with enzymes of the pathway, a data-driven iterative modeling approach was applied. In this cycle, an abstract model of the β-ketoadipate (ortho-cleavage) pathway was refined until the measured data were described mathematically. The most appropriate model resulting from this iterative model-building cycle is outlined in Fig. 1, with the corresponding mass-balance and kinetic rate expressions. In total, the model comprises 27 timedependent variables, of which 16 are states expressed as ordinary differential equations (Table 1), and 11 are reactions involved in either the transformation or transport of metabolites (Table 2). Multiple assumptions were essential to reduce the model complexity, especially in the case of insufficient in vitro or in vivo information regarding the enzymatic reactions. Some cardinal assumptions that were made to model the steady-state and transient behavior of the system are as follows: (a) the biomass, growth rate, and cell volume do not change during the transience after perturbation, (b) the mixing in the reactor is homogenous, and (c) a constant external pH of 7 and a membrane potential of À 50 mV are maintained. Although the pulsed benzoate was consumed completely in an overall time span of $25 min, the pulse response methodology applied for model validation was restricted to 12 min after perturbation. Within this time window, we assumed that the changes were primarily caused by metabolic regulation and by changes in biosynthetic reactions (e.g., protein biosynthesis was regarded as remaining in a "frozen" state) . Consequently, the maximal enzyme capacities were assumed to remain at a steady-state value (for calculation ofr max values, see Section 3.4). Apart from these intrinsic assumptions, the biochemical assumptions made to describe the kinetic properties of the enzymes are detailed in Section 3.3.
The time-series metabolome data generated from benzoatepulse experiments were integrated into the model to identify the kinetic parameters using a stepwise internalization approach that was proposed previously . For this purpose, the measured concentrations of intracellular benzoate cis-diol ( ) C DIOL in and intracellular catechol ( C CAT in ) and the extracellular concentrations of benzoate ( C BEN ex ), catechol ( C CAT ex ), cis,cis-muconate ( C MUC ex ) and β-ketoadipate (C KA ex ) were included in the model identification procedure. Additionally, for metabolites with high extracellular concentrations, mass-balance equations were derived to introduce their total concentrations (C ; C MUC tot KA tot ) into the model structure. The measurements were then numerically integrated with simulated concentrations from the mass-balance equations, and through minimization of the squared error, the unknown parameters were estimated (Table 3). For the model that resulted from this identification procedure, the predictions were in satisfactory agreement with the experimental observations. The identified model was able to predict both the steady-state and transient conditions introduced by a benzoate pulse (Fig. 3). To further validate the predictions of the constructed kinetic model, oxygen balance equations for the liquid and gas phases were included. The model fit the observed trend in the measured values of oxygen in the gas and liquid phases without additional parameter estimation (Fig. 4).

Thermodynamic feasibilities of benzoate metabolism reactions
Due to the limited information regarding the regulatory interactions of the benzoate metabolic network, we examined the feasible ranges of benzoate metabolic network operation via thermodynamic analysis. According to the second law of thermodynamics, a reaction proceeds in the direction of negative ∆ G r . NET analysis (Zamboni et al., 2008) was performed to elucidate the active regulatory sites of benzoate metabolism. NET analysis couples metabolite concentrations to an operating metabolic network and checks for thermodynamic feasibility based on the second law of thermodynamics (i.e., to allow for negative Gibbs energy for each preset flux). The steady-state intracellular concentrations of benzoate cis-diol and catechol and the extracellular concentrations of catechol and β-ketoadipate were used for NET analysis. From a flux balance analysis tool that is integrated into the anNET software, an operating benzoate metabolic network Fig. 4. Comparison between measurements (filled circles) and model predictions (continuous lines) for oxygen in the liquid (top) and gas phases (bottom) during the benzoate pulse (at time zero) response experiment in a steady-state benzoatelimited culture. The subscript 'RV' on the y-axis represents the reactor volume. The solid lines are the results predicted using the dynamic model without additional parameter adaptation. Thus, these lines serve to validate the dynamic model. with 154 reactions and 183 metabolites was deduced from the P. putida genome-scale metabolic model iJN746 (Nogales et al., 2008). Thermodynamic consistency was demonstrated for the operation of the benzoate metabolic network with the given metabolite measurements, and the concentration ranges of the unmeasured metabolites were calculated for the feasible ranges of the estimated Gibbs energy of reactions. Of the 154 reactions of the benzoate metabolic network, the reactions of the β-ketoadipate (ortho-cleavage) pathway, the tricarboxylic acid cycle, gluconeogenesis, and the non-oxidative pentose phosphate pathway were analyzed in detail (Fig. 5).
Further, in order to extend the link between thermodynamic analysis and parameter sensitivities, additional analysis has been performed by calculating the flux control coefficients with the software SBtoolbox2. The matrix of flux control coefficients is shown in (Fig. 6). The flux control coefficients can be linked to the free energy via the elasticities. Westerhoff et al. (Westerhoff et al., 1987)  i.e. the reaction catalyzed by benzoate cis-diol dehydrogenase (K mDIOLDH , K mCATDH , K iDIOLDH ), catechol 1,2-dioxygenase (K mCDO , n cO2 ), and β-ketoadipate succinyl-CoA transferase (K mKACOAT ) were found to exert significant flux control over the pathway and the kinetic parameters of the reactions with low ∆ G r have no flux control over the pathway.

Discussion
We investigated the well-known but poorly understood in vivo regulatory interactions of the β-ketoadipate (ortho-cleavage) pathway in P. putida KT2440 using an iterative model building process. Critical assumptions were necessary to describe the in vivo behavior of metabolites with appropriate reaction kinetics. In addition to metabolic regulation, there must be an additional regulation phenomenon at the enzyme level; this phenomenon was observed during the entire course of our pulse response experiments. The changes in metabolic regulation is known to occur on the time scale of seconds . The hypothesis that an additional regulation phenomenon is superimposed on metabolic regulation was supported by the increasing degradation rate of benzoate with time (Fig. 7A). This observation is the main reason for restricting the metabolic response to 12 min after perturbation and not to a longer scale (12-23 min). The importance of additional regulation at the enzyme level was further illustrated in an experiment with sequential pulse experiments ( Fig. 7B and C), in which the uptake rate of benzoate increased as the number of generations increased in a continuous culture kept at a constant dilution rate. Although the linear increase in the substrate uptake rate is a characteristic feature of continuous cultivation, whether the entire pathway or only the benzoate transport enzyme is optimized for carbon uptake is unknown and must be investigated thoroughly. A previous study showed that when the yeast S. cerevisiae is confronted with multiple glucose stimuli, the specific uptake rate of glucose is increased by a factor of seven compared with its steady-state value; this change accompanies altered expression of glucose transport proteins (Buziol et al., 2008). Moreover, prolonged glucose limitation of S. cerevisiae in a chemostat causes multiple steady states with different metabolite and enzyme levels (Wu et al., 2006). The observed deviations between Fig. 6. Flux control coefficients. The colors in row i and column j indicate the control that the parameter j exerts over the reaction rate i. The colors denote positive and negative control, respectively. A small increase in the Km value of catechol (i.e. K mCDO ) or the hill coefficient (n cO2 ) exerts a negative control over all the fluxes of the pathway, and a positive control over the transport flux of benzoate cis-diol ( r DIOLT ). With respect to the enzyme benzoate cis-diol dehydrogenase, a small increase in the values of kinetic parameters K mDIOLDH , and K mCATDH exerts a negative control over the transport flux of benzoate cis-diol (r DIOLT ) and a positive control over the transport flux of βketoadipate ( r KAT ). A small increase in the value of K iDIOLDH exerts a negative control over all the fluxes of the pathway and a positive control over the transport flux of benzoate cis-diol (r DIOLT ). With respect to the enzyme β-ketoadipate succinyl-CoA transferase, a small increase in the Km value of this enzyme (K mKACOAT ) exerts a significant negative control over the transport flux of β-ketoadipate ( r KAT ). Also, a small increase in the kinetic parameter of cis-cis muconate (K mMUCT ) and β-ketoadipate (K mKAT ) transporters were found to exert a positive control over the transport rate of cis-cis muconate (r MUCT ) and β-ketoadipate (r KAT ). the model predictions and the measurements (Figs. 3 and 4) also indicate the presence of an additional regulation phenomenon in P. putida that cannot be compensated by parameter fitting with the existing model structure.
Major assumptions were made with respect to the feedback inhibition of benzoate transport and benzoate cis-diol dehydrogenase. Although no evidence to support a molecular basis for this phenomenon was reported in Pseudomonas for the former case of benzoate transport inhibition by benzoate cis-idol, a study with Corynebacterium glutamicum has shown that hydroxyl-substituted benzoates inhibit benzoate uptake significantly (Wang et al., 2011). In the latter case, feedback inhibition of benzoate cisdiol dehydrogenase by catechol was used purely to describe the in vivo dynamics of benzoate cis-diol and catechol concentrations. In addition, the reason for the low in vivo concentration of catechol compared with benzoate cis-diol might be the existence of two dedicated catechol dioxygenases, of which CatA2 was recently described as a safety valve that protects against toxic concentrations of catechol (Jimenez et al., 2014). Consequently, all of these factors might result in significant control of the total pathway flux by the benzoate cis-diol dehydrogenase.
These five reactions were classified as non-equilibrium reactions because at large À ∆ G r , the thermodynamic driving force is significant, and the rate of the reaction is regulated predominantly by changes in the enzyme capacity (for example, via transcription or the binding of effectors). These indications for sensitivity to additional regulation phenomena are partially manifested in the corresponding kinetic expressions. In our benzoate-pulse experiments, the substrates of these reactions were also found to be extracellular; this observation was true to a greater extent for some metabolites, such as β-ketoadipate, benzoate cis-diol, and cis,cis-muconate. The excretion or the overflow of the individual metabolites is most likely a result of the adaptation of the culture to the steady-state conditions of the continuous culture. Some of these excreted metabolites are extremely toxic (e.g., catechol), and other metabolites may be excreted due to simple metabolic overflow. The scenario can be changed by using different dilution rates, which may result in different concentrations of the enzymes. However, these phenomena were not studied in this project.
Bacterial mineralization of aromatics is well known to be critical e.g. on the level of benzene ring cleavage (Dagley, 1971;Evans, 1963). We initially hypothesized that catechol accumulation in shake flasks experiments are regulated by ring cleavage enzymes. The model and especially the quantitative description of all the rates allowed us to verify this hypothesis. In addition, catechol formation and especially benzoate uptake were identified as key and previously not recognized regulation points. For benzoate uptake, this new hypothesis was experimentally verified by sequential pulses. Amazingly cells reacted via aging/adaptation to the bottleneck benzoate uptake by increasing the benzoate uptake rate ( Fig. 7B and C). This now can be taken as starting point again for defining new objectives, e.g. substrate uptake, product formation, etc. to run the systems biology iterative process (experiment-model-hypothesis-experiment…) a second time. In addition, model building is a purpose-driven-process. In other words: first define the purpose before starting to build the model. The purpose for the presented dynamic model is the application of re-designing the enzymatic makeup for the process/product of interest. In our specific case, objectives for follow up biotechnology research might be rate improvements for benzoate uptake and catechol formation and cleavagethe major flux bottlenecks identified by in vivo dynamic modeling (Fig. 1) on whatever level appropriate (metabolic engineering, reaction engineering, process engineering).
In summary, a data-driven modeling approach increased our understanding of the complex interactions of the β-ketoadipate (ortho-cleavage) pathway. The proposed model can capture the stationary and dynamic behavior of intracellular and extracellular metabolites. The systematic integration of quantitative metabolite data with thermodynamic information was pivotal for demonstrating that only a few reactions of benzoate metabolism are candidates for active regulation (e.g., transcriptional and/or allosteric regulation), whereas the rates of most reactions depend only on substrate concentrations that are close to equilibrium. Although the critical assumptions that were made regarding the regulation of the benzoate transport rate and benzoate cis-diol dehydrogenase activity are supported by in vivo dynamic modeling, they do invite future confirmation by in vitro investigations of these proteins. These findings highlight the necessity of quantifying metabolic rates of well known metabolic pathways in vivo.
The biochemistry of pathway enzymes in the cytoplasm has to be described both theoretically in (dynamic) models as well as verified experimentally, e.g. in chemostat experiments. In this respect, direct translations of genomic sequence data/ annotations into metabolic flux networks are necessary, but not sufficient.