An energetic reformulation of kinetic rate laws enables scalable parameter estimation for biochemical networks

The technology for building functionally complete or ‘whole-cell’ biological simulations is rapidly developing. However, the predictive capabilities of these simulations are hindered by the availability of parameter values, which are often difficult or even impossible to obtain experimentally and must therefore be estimated. Using E. coli ’s glycolytic network as a model system, we describe and apply a new method which can estimate the values of all the system’s 102 parameters – fit to observations from studies of proteomics, metabolomics, enzyme kinetics and chemical energetics – and find that the resulting metabolic models are not only well-fit, but also dynamically stable. An analysis of how well parameter values in the network were determined by the training data revealed that over 80% of the parameter values were not well-specified. Moreover, the distribution of well-determined values was biased to a specific part of the network and against certain types of experimental data. Our results also suggest that perturbing the functional, energetic space of parameters (rather than traditional metabolic parameters) is a superior strategy for exploring the space of biological dynamics. The estimated parameter values matched both training data and previously withheld validation data within an order of magnitude for over 85% of the data points; notably, the area of greatest frustration in the network was also the most fully determined. Finally, our estimation method showed that fidelity to physiological observations such as network response time is enforced at the cost of fit to molecular parameter values. In summary, our reformulation enables estimation of accurate, biologically relevant parameters, generates insight into the biology of the simulated network, and appears generalizable to any biochemical network – potentially including whole-cell models. © 2018 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
In 2012, the construction and output of a 'whole-cell' model was reported for the simplest known culturable bacterium, Mycoplasma genitalium ( Karr et al., 2012 ). In brief, a whole-cell model is an integrated simulation of several 'processes' or dynamical submodels, which attempts to account for the known function of every gene in a given organism over at least one cell cycle. The separation of the model in this way -for example, into one submodel that represents transcription and another that represents carbon and energy metabolism -enables a process-specific selection of modeling paradigms in order to best capture the available data relevant to that process as well as the behavior which is to be mod-eled. Subsequent analysis of the M. genitalium model showed that it was able to predict the results of previously unperformed experiments, both at the level of cell physiology as well as enzyme kinetic parameters . A number of follow-up studies have built on this foundation, providing new applications ( Purcell et al., 2013 ), databases ( Karr et al., 2014;, and visualizations ( Lee et al., 2013 ).
As large-scale in silico proxies for real cells, whole-cell models have the potential to serve as platforms for integrating diverse sources of biological data and generating new, testable hypotheses consequent of observed model behavior ( Carrera and Covert, 2015 ). However, the assembly of a complete whole-cell model is a daunting task ( Macklin et al., 2014 ). Throughout the iterative refinement of submodels, progress is restricted by the lack of specific quantitative data needed to parameterize a known dynamical behavior. Further, quantitative data may be inconsistent; this is particularly a concern when data is gathered from different stud- ies, on different organisms, or in non-native environments (e.g. in vitro rather than in vivo ). An excellent review of the particular challenges posed by whole-cell model parameterization was published last year ( Babtie and Stumpf, 2017 ). These challenges limit the scope of whole-cell models, and subsequently, their predictive power. To support the continued refinement of whole-cell models, it is clear that we need to adopt and develop strategies for wholecell parameter estimation.
Towards this end, we sought to develop parameter estimation techniques applicable to whole-cell like models, starting first with the intention of expanding the metabolic submodel. To date, the whole-cell metabolic submodels have utilized variants ( Birch et al., 2014 ) of flux-balance analysis (FBA), a constraint-based optimization technique that strictly enforces stoichiometric limitations but does not directly incorporate kinetic or energetic information ( Orth et al., 2010 ). FBA is advantageous in that it can simultaneously capture the qualitative functional role of thousands of genes with relatively few model parameters. However, the top-down, optimization-based nature of FBA limits its generalizability, as the constraints must be manually tuned for each modeled condition. Alternatively, a traditional metabolic model composed of kinetic rate laws in a system of ordinary differential equations (ODEs) would allow for bottom-up, mechanistic predictions of metabolic activity. As a model system for the development of parameter estimation techniques, we selected the glycolytic pathway of E. coli both for its detailed characterization in the literature and its immediate relevance to our ongoing work on the E. coli whole-cell model. ODE models of biochemical pathways require exquisite, complete parameterization of reaction activity; our model, composed of kinetic rate laws assembled using a procedure analogous to convenience kinetics ( Liebermeister and Klipp, 2006 ), is no exception. To inform our parameter values, we utilized observations from studies of proteomics, metabolomics, chemical energetics and enzyme kinetics. However, as previously discussed, the available data is neither complete nor is it wholly applicable to in vivo models. Therefore obtaining a complete, consistent set of parameter values requires us to formulate and solve an explicit parameter estimation problem, with the goal of matching some desired dynamical behavior while simultaneously minimizing disagreement with the available observations.
These parameter estimation problems typically take the form of multi-modal optimization problems (i.e. problems with multiple local 'solutions'); consequently they require global optimization algorithms. Global optimization algorithms can be quite computationally expensive (i.e. running for hours on hundreds or thousands of processors), and furthermore there is no guarantee for convergence ( Moles et al., 2003 ). Despite these challenges, a number of successful large-scale metabolic parameter estimations have been performed. Most analogous to the work presented here are several detailed kinetic models of E. coli central carbon metabolism ( Jahan et al., 2016;Kurata and Sugimoto, 2018;Millard et al., 2017 ), optimized using traditional approaches like genetic algorithms. Other approaches (such as the ORACLE Chakrabarti et al., 2013 andk-OptForce Chowdhury et al., 2014 platforms) utilize FBA to generate flux profiles in concert with kinetic model parameterization. Additionally, an 'ensemble modeling' approach ( Khodayari et al., 2014 ) eschews typical Michaelis-Menten-like kinetic rate laws in favor of parameterizing elementary reactions. Beyond parameterizing these kinetic models, these large-scale parameter estimation effort s have yielded insights into the systems-level control of metabolic activity. However, indicated throughout these works is a consistent need for powerful, efficient estimation of parameters. Parameter estimation is itself an active area of research, including both statistical (e.g. likelihood maximization, Bayesian inference ( Hug et al., 2013;Liepe et al., 2015 )) and optimization-based approaches ( Fröhlich et al., 2017;Jagiella et al., 2017;Penas et al., 2017 ). Because our eventual goal involves whole-cell scale estimation, which is currently beyond the capacity of statistical approaches ( Babtie and Stumpf, 2017 ), we focused on optimization for the work described here.
In the development of our own parameter estimation techniques, designed with needs of whole-cell-like models in mind, it was found that a reinterpretation of kinetic rate laws (our underlying dynamic models) in terms of energetics facilitated several downstream procedures and provided a new conceptual basis for interpreting the space of metabolic dynamics. The details of this reinterpretation and its consequences are the subject of this study.

Results
In Section 2.1 we introduce our energetic reformulation of standard kinetic rate laws, used throughout the study. In Section 2.2 this new formulation is used to connect observations to parameter values, elucidating the limited extent of existing literature in parameterizing a whole-cell-like kinetic model. We further apply this formulation to derive a new approach for exploring the parameter space, dubbed 'parsimonious' perturbation, in Section 2.3 . The resulting model parameters and their fit to training data are interpreted in Section 2.4 , with the dynamics of the parameterized models following in Section 2.5 . We conclude with a discussion about the ramifications of this new formulation as well as the biological consequences of our optimized parameter values. Short explanations of our approaches are provided in the Appendix. Full mathematical derivations, explicit methods, secondary results and a complete record of all curated data are in the provided Supplementary Materials.

Standard kinetic rate laws can be reformulated to an energetic basis
For this study we developed an energetic reformulation of Michaelis-Menten-like kinetic rate laws, which greatly facilitates downstream parameter estimation procedure and unifies several energetic ideas (e.g. transition energies and binding energies) into one Boltzmann distribution-like equation. The most immediate utility of this reformulation is its decomposition of several coupled parameters (e.g. forward and reverse catalytic rate constants, equilibrium constants) into uncoupled (i.e. independent) parameters (e.g. standard energies of formation). We also discovered several additional benefits to this reformulation; a handful of these will be described in subsequent sections, with many others expanded in the Supplementary Materials.
Our reformulation begins with a simple metabolic model expressed as a system of ordinary of differential equations (ODEs). These equations describe the time evolution of intracellular metabolite concentrations ˙ c with respect to the stoichiometries of modeled reactions S and kinetic rate laws for each reaction v (c) (which are themselves functions of the metabolite concentrations). Additionally, we account for the dilution of metabolites (consequent of cell growth) at the specific growth rate μ. Together this system of ODEs can be compactly represented as The kinetic rate laws v ( c ) are assembled using a procedure that yields equations comparable to convenience kinetics, an approach developed by Liebermeister et al. just over a decade ago ( Liebermeister and Klipp, 2006 ). Briefly, every kinetic rate law is designed to be 'complete' in the sense that it is fully reversible (subject to the limitations imposed by chemical energetics) and exhibits proper asymptotic behavior (i.e. saturating to some maximum rate at high substrate concentrations, and likewise diminishing to a rate of zero at low substrate concentrations). We use this procedure to assemble an ODE model of 102 parameters for eleven kinetic rate laws and eighteen metabolites. The full kinetic rate laws for the model system can be found in the Supplementary Materials.
For the sake of demonstration, we consider the kinetic rate law associated with the simple case of a reversible reaction This is a standard formulation of Michaelis-Menten-like kinetics for a unimolecular, reversible reaction (i.e. an isomerization). This formulation is appealing in that it is naturally expressed in terms of easily interpreted parameters (e.g. concentrations, catalytic rate constants). In contrast, our reformulation expresses these kinetic rate laws in terms of energetic parameters. Under this transformation, developed in Appendix A.1 , Eq. (2) becomes where E( G ) = exp − G/RT , k is an arbitrary constant for correcting units, and are the differences in energy levels ( G terms) as linear combinations of other energy or energy-like parameters. The relationships between these basic and derived parameters are depicted schematically in Fig. 1 . Our reformulated parameters can be interpreted in terms of existing kinetic and energetic concepts; for example, G MAK,f is the transition energy of the forward reaction in the limit that the rate of reaction is well approximated by mass-action kinetics (i.e. when substrate concentrations are relatively low), and G b A is the energy associated with the propensity for molecule A to bind the enzyme. Likewise the new basic parameters are related to well understood terms; G • A is the standard Gibbs energy of formation for molecule A, and G c A is the Gibbs energy contribution associated with concentration (i.e. G c A = RT ln c A ). (More thorough descriptions of these terms and their relationships with other energetic or kinetic parameters are provided in the Supplementary Materials.) In the following sections we will describe how this reformulation facilitates both basic analysis of our training data as well as the development of our novel perturbation approach for parameter estimation.

The energetic reformulation readily links data availability to model parameters
The parameter reformulation presented in Section 2.1 allows us to ascertain the degree to which our model parameter values can be determined from training data, even when observed parameters and model parameters do not have a one-to-one relationship (e.g. in the case of chemical equilibrium constants and standard energies of formation). The data used in this study is described in Appendix A.2 , while the description of our 'determinacy analysis' procedure is given in Appendix A.3 . In brief, the training data were transformed into the energetic parameter space (as a vector b ) and associated with the basic parameter values (as vector x ) via the linear expression b = F x ; the matrix F can then be analyzed to determine where a parameter value is undetermined (the training data provides no information to help determine the parameter value) fully determined (the value can be completely determined from the training data), or partially determined (the value is informed by the training data, but requires further information to be fixed). The determinacy of each parameter is depicted on our network model of glycolysis in Fig. 2 .
In applying this procedure to our model system, we find that, of the 102 model parameters, 19 are fully determined, 32 are partially determined (i.e. only determined relative to other parameters), and 51 are wholly undetermined. Thus, over 80% of the parameters in our model are underdetermined.
Beyond the total determinacy of parameter values, given the curated literature, we also observe biases in the availability of certain types of data. First, we observe that part of the network, centered around the fructobisphosphate aldolase (fba) reaction in upper glycolysis, has been investigated more extensively than the rest. Additionally, the majority of underdetermined parameter values (35/51 or 69%) correspond to effective binding energies (i.e. Michaelis-Menten saturation constants), particularly for substrates of the reactions opposite the canonical (usually glycolytic) direction of reaction.
We also see that both enzyme concentrations and standard energies of formation are never fully determined. This is readily explained in the context of the training data. In the case of standard energies of formation, these parameter values are only determined relative to one another, due to the fact that our most relevant training data is chemical equilibria constants (each of which is a composition of two or more standard energies of formation in their corresponding stoichiometric ratios). The underdeterminacy of the enzyme concentrations arises from the fact that proteomics data were expressed as relative abundances, eliminating the need to scale this data but somewhat reducing its information content.  Fig. 1 B). Each symbol is styled according to its determinacy , i.e. the extent to which its value can be inferred from the training data (parameter value observations).
While not all independent parameters need to be established to concretely define model dynamics, the extent of non-fully determined parameters, as well as the non-uniformity of where determinacy occurs throughout the network, demonstrates a serious need for a parameter estimation system that can estimate a large number of unknown parameter values simultaneously. In the next section we will perform such a parameter estimation, facilitated by our energetic parameter reformulation.

A novel perturbation approach leads to more successful parameter estimation
Next, we recognized that the energetic reformulation of the kinetic rate laws enabled the development of a novel perturbation approach that leads to successful estimation of parameter values, even when applied to a rudimentary optimization algorithm. We began with a 'naive' approach in which the basic model parameters were perturbed individually. This approach struggled to converge to a dynamically viable, well fit set of parameter values, as confirmed by inspecting objective values as well as performing linear stability analysis ( Strogatz, 1994 ) (LSA) on the parameterized models (see Appendix A.6 ) ( Fig. 3 ). This result was illustrative of the technical challenges posed by these parameter estimation problems. Several strategies have been employed in the context of biological networks, with varying degrees of success ( Moles et al., 2003 ). In practice, strategies such as genetic algorithms ( Jahan et al., 2016 ) and particle swarm optimization ( Millard et al., 2017 ) have been successfully applied to estimate parameter on models of large biochemical networks, albeit at significant computational costs.
These results led us to suspect that a perturbation strategy which was sensitive to the context of the parameters (and therefore sensitive to the structure of the modeled biological system) could more effectively search the parameter space. For example, consider a parameter such as the concentration of ATP, a metabolite utilized in several reactions. Somewhat near an optimal set of parameter values, it would be difficult to perturb the concentration of ATP without upsetting the dynamics in many other parts of the network, presumably causing the optimization to become trapped in a 'corner' of the high-dimensional parameter space.
To address this concern we developed an approach that allowed us to perturb the parameters in a fashion that was minimally disruptive to the dynamical behavior of the system as a wholehence, a 'parsimonious' perturbation in terms of system behavior. Parsimonious perturbations can be achieved by perturbing derived model parameters related to dynamical behavior (e.g. the differences in energy levels corresponding to a maximum rate of reaction v max or a saturation ratio c/K M ), and selecting a search direction that minimally disrupts the other derived model parameters. The mathematical procedure for obtaining a parsimonious perturbation vector is given in Appendix A.4 .
For either perturbation approach, our parameter estimation problem is formalized as the simultaneous optimization of two objective terms; first, to minimize the disagreement with reported parameter values, expressed as the 'misfit cost' function f ( x ); and second, to ensure that the simulations achieve a steady-state, which we express as the 'non-steady-state cost' function g ( x ). Formally we write this as the optimization problem where w is some weight that we increase until g ( x ) is sufficiently small. (The complete expression of these terms is provided in Appendix A.5 .) To solve this optimization problem, we employ an algorithm (described extensively in the Supplementary Materials) in which we perturb our current best parameter values and search for an improved (smaller) objective value. The resulting values of these objective terms are shown in Fig. 3 for 300 attempts using each perturbation approach. Fig. 3 shows that, in contrast with the 'naive' approach, the majority of attempts using the 'parsimonious' perturbation approach successfully minimize the non-steady-state cost (234 or 78%). These attempts also exhibit misfit costs that are at the lowest part of the range. Our results suggest that perturbing the functional, energetic space of parameters (rather than traditional The values of the two major objective components following three-hundred (300) optimization attempts, using the 'naive' and 'parsimonious' perturbation approaches. The dynamic 'viability' of each parameter set was determined via simulation of the parameterized models and linear stability analysis. 234 (78%) of the attempts using the 'parsimonious' approach were found to be viable, compared to just 2 (0.67%) of those using the 'naive' approach.
metabolic parameters) is a superior strategy for exploring the space of biological dynamics.

Estimated parameter value distributions compare well with training and validation data
Our next objective was to evaluate the quality of the dynamically viable sets of estimated parameter values, using three metrics: (1) fit to the training data, (2) fit to a set of withheld validation (test) data, and (3) consistency of the output simulation dynamics with physiological observations. Fig. 4 B compares the estimated parameter values against the 144 training data values (corresponding to 45 basic or derived model parameters) utilized in the expression of the misfit cost function. We find that the median estimated parameter values are, as desired, generally well fit to their training (40 of 45 or 88% falling within ten-fold of the median observed value), and moreover that the distributions of the parameter values are typically quite narrow.
Beyond the degree of fit to our training data, we also sought to evaluate whether our parameter estimation system could make reasonable predictions of withheld parameter values (i.e. test or validation data). Fig. 4 C shows the distributions of predicted parameter values against catalytic rate constants estimated from specific activity data. Here we find that our predictions are 'well fit' (again, within ten-fold) with the observed values for five of the six reactions, where the one poor fit is for the phosphoglycerate mutase (gpm) reaction.
The main point of frustration in these comparisons appears to be centered around the triosephosphate isomerase (tpi) reaction, a topologically unusual reaction for this network in that it interconverts between the two upstream products of the fructobisphosphate aldolase (fba) reaction to feed into the glyceraldehyde phosphate dehydrogenase (gap) reaction. Parameter values related to these three reactions are poorly fit to the training data, suggesting one or more conflicts between the experimental data and the desired behavior of the dynamical model. In particular we find that the chemical equilibria constants for the fba and tpi reactions, along with the concentration of the fba enzyme and Michaelis-Menten saturation constant for glyceraldehyde phosphate in the gap reaction are all shifted in a fashion that would facilitate greater flux through glycolysis. Interestingly, this portion of the network is also the most fully determined with regard to model parameters.

Identifying a trade-off between conformity to physiological observations and fit to molecular-level measurements
We next considered our third metric listed above: the consistency of the output simulation dynamics to physiological observations. We first verified that all sets of estimated parameter values with sufficiently low non-steady-state costs (very small values of g ( x ), as shown in Fig. 3 ) were also dynamically stable, reaching a steady-state over time (top row of Fig. 5 ). However, this characterization also revealed that the recovery times were very slow, on the order of hours. This was far longer than the anticipated perturbation recovery times for real cells, which other studies suggested should be on the order of seconds ( Schaub and Reuss, 2008 ).
By performing LSA on the original model equations ( Eq. 1 ), under the assumption that all metabolite concentrations were in excess of their corresponding Michaelis-Menten saturation constants ( K M values), we found that all characteristic recovery times (CRTs) converged to τ = 1 /μ ≈ 1 . 4 h , consistent with the re- covery times observed by simulating the ODEs. This suggested that the cause of the slow recovery times was high enzyme saturation.
To test this hypothesis, we augmented the original parameter estimation problem to include terms in the misfit cost function f ( x ) that penalized for excess enzyme saturation. While the misfit cost function was originally designed to penalize for order-ofmagnitude differences between training data and parameter values, our energetic reformulation of model equations allows us to also penalize for any combination of parameters values against some other expected value. We chose to penalize for metabolite concentrations c in excess of their corresponding Michaelis-Menten saturation constants K M (i.e. when c/K M > 1 ). The parameter estimation procedure was repeated on this augmented problem, using a range of weights on the new terms.
As the weighting on these terms was increased, we found that the rates of recovery following a perturbation increased dramatically (as shown in the remaining panels in Fig. 5 ). This is partic-ularly apparent when the weight on the terms is greater than the weight on the training data ( 10 +1 and higher).
We also found that reduced enzyme saturation (higher weights on the saturation penalty) led to other effects on our parameter estimation, summarized in Table 1 . First, we noted that the overall determinacy of the network was increased considerably by introducing our penalty terms for enzyme saturation. In particular, whereas 50% of the parameters were undetermined previously, our new strategy reduced this percentage to 9%. The inclusion of enzyme saturation constraints also led to more realistic simulations in terms of pathway utilization, as the canonically glycolytic reaction pyruvate kinase (pyk) was favored over the canonically gluconeogenic reaction, phosphoenolpyruvate synthetase (pps). In contrast, the low-and no-weight parameter estimates tended to favor the pps reaction.
The benefits of higher network determinancy and more realistic simulations with respect to response time and pathway utilization were tempered by an increase in the misfit cost. Examining the parameter misfit in more detail, we found that seven param-  eter values were impacted, again centered about upper glycolysis ( Fig. 6 A). The degree of fit to the catalytic rate constants for pyk and pps in our validation data set were also negatively impacted ( Fig. 6 B). This was unsurprising, given that the training data had originally favored the pps reaction over pyk as described above.

Discussion
In summary, we found that a reformulation of kinetic rate laws in terms of energetic parameters was necessary for developing parameter estimation techniques applicable to whole-cell modeling. Such a reformulation was partially considered in the development of convenience kinetics ( Liebermeister and Klipp, 2006 ); our approach builds on those concepts by also reformulating metabolite and enzyme concentrations as energy or energy-like quantities. While care should be taken to not overly interpret these quantities as true energies (e.g. G b should be most properly thought of as an effective binding energy), this energetic reformulation nonetheless provides a theoretical basis for exploring the space of biological dynamics in the derived, functional coordinate space of differences between energy levels.
An interesting insight that arose from the energetic reformulation came as we analyzed the determinacy of the network, using a technique that can link parameter value observations to model parameters, even when the relationship between the two is complex. In the context of this study, we found a considerable amount of underdetermined parameter values, consequent of a lack of applicable experimental observations. In some cases the observations may have been made at some time but were excluded by our somewhat selective curation process (e.g. our measurements were required to come from studies of E. coli ; see Appendix A.2 for more requirements). For other parameters, no measurements have been made to date. In any event, we found that the determinacy was not only biased by experimental method but also by location in the network. Thus, our determinacy analysis should help in prioritizing future measurement technologies to be developed and experiments to be performed. However, considering that our system is one of the most well-established pathways in perhaps the world's most characterized organism, this finding underscores the need for large-scale parameter estimation techniques to better understand and model biological systems.
Most critically, we discovered that our energetic reformulation enabled us to develop a new perturbation approach, which enabled a rudimentary optimization algorithm to reliably estimate well fit, dynamically viable model parameters. In particular, we found that perturbing the functional, energetic space of parameters (rather than traditional metabolic parameters) is a superior strategy for exploring the space of biological dynamics, leading to sets of parameter values that fit the training and validation data well and were also dynamically viable.
We were surprised to observe a slow response time in our original simulations (i.e., based on the first parameter value estimates). We addressed this issue by motivating our parameter estimation system to find solutions with reduced enzyme saturation. While we may naturally expect enzymes to be heavily saturated, thereby achieving the highest possible rate of reaction per enzyme, our results suggest that this comes at the cost of reduced adaptability to fluctuations in metabolite concentrations. This further implies a competing evolutionary pressure between enzyme efficiency and enzyme sensitivity, i.e. a trade-off between efficient growth under stable conditions and flexible growth under variable conditions. Whether this trade-off is a generalizable property of (unregulated) metabolic pathways remains to be seen.
Our reduced enzyme saturation approach led to more realistic system-level behaviors, including faster recovery following a perturbation as well as utilization of the more canonical pathway for pyruvate production. However, these improvements came at the cost of worse fits to training and test data. These results can be explained in three ways: first, that kinetics data for glycolysis -typically generated in vitro -is not consistent with expected in vivo behavior, due to the influence of ions, pH, and molecular crowding in the cytoplasm that could lead to different effective parameter values ( Bennett et al., 2009 ). Alternatively, we could speculate that our expectations about pathway utilization and perturbation recovery times are incorrect. A third explanation, which requires discounting neither the parameter value observations nor our behavioral expectations, is that our saturation penalty partially captures the effects of kinetic (e.g. allosteric) regulation (effects which are not explicitly included in this model). Like the increased sensitivity conferred by the reduced levels of saturation, allosteric regulation would allow modulation of reaction rates subject to perturbations in the system (however these would likely be global signals, rather than the local inputs and outputs of a given reaction). Glycolysis is in fact known to be allosterically regulated ( Frankel, 1996 ), lending credence to this possibility.
Finally, our parameter estimation also proved capable of identifying and reconciling points of frustration in the training data. Of particular note was our observation that the chemical equilibria constants associated with the fructobisphosphate aldolase (fba) and triosephosphate isomerase (tpi) reactions had to be shifted to favor the forward (glycolytic) direction of reaction to achieve sufficient flux through the pathway. This result is parallel to our findings regarding the catalytic rate constants for pyk and pps, in which adding the constraint on enzyme saturation unexpectedly led to simulations that used the canonical glycolytic pathway -but at the cost of agreement with some validation data. We expect that further targeted measurements may be able to resolve these frustration points.
We believe that the approach described here should generalize to any biochemical network (as it is informed by the equations, which are in turn informed by the network structure) as well as to other optimization algorithms. In fact, we see no reason why this technique could not generalize to domains beyond biology, so long as systems of equations can be expressed as sums of polynomials, rational functions, or Boltzmann distributions.
Ultimately, a complete whole-cell model should simulate the functional role of every gene and gene-product. This poses a tremendous parameter estimation challenge, on a scale far beyond the simple pathway we have modeled in this study. Even neglecting regulatory details, we estimate that a complete kinetic model of metabolism, including hundreds of metabolites and reactions, would easily require thousands of parameters. Future work will be required to conclusively determine whether the methods described here can truly be expanded to that scale, and even beyond metabolism to other biological submodels. If so, we anticipate that one of the most challenging bottlenecks in whole-cell modeling can be largely resolved.

Data declaration
All code and data used in this study are available at https://github.com/CovertLab/parest.
where G b A is the effective Gibbs energy for the binding of metabolite A to the enzyme. (For expressive convenience, and without loss of generality, we have selected a standard reference energy of 0 kcal mol −1 for a concentration of 1M.) Substituting these transformations into Eq. (2) , the kinetic rate law becomes where we define E(x ) = exp −x/RT for the sake of brevity. We can further rearrange this expression to make the forward and reverse reaction rate terms in the numerator symmetric: Finally, we redefine the forward catalytic rate constant as where G ‡ is defined as the energy of the transition state, and k is simply a proportionality constant that could be estimated from the Eyring equation ( Eyring, 1935 ). Together these transformations give us Eq. (3) . Additional relations, as well as a glossary of symbols, can be found in the Supplementary Materials.

A2. Data curation
The pathway reaction substrates, products, enzymes and stoichiometry were obtained from EcoCyc ( Keseler et al., 2017 ). A complete description of the curated training data (sources and values) are provided Supplemental Materials. In brief: Intracellular metabolite concentrations were mostly taken from Bennett et al. (2009) . Some data were excluded; for example, the hexose-6-phosphate concentration, which would include both glucose-and fructose-6-phosphate, was not used because of its ambiguity. The proton concentration was computed assuming a pH of 7.4 ( Wilks and Slonczewski, 2007 ). The water concentration was computed assuming a water content of 70% ( Neidhardt et al., 1990 ). Altogether, concentration training data for eleven of eighteen metabolites were used. In addition to the concentration data, we also utilized 61 chemical equilibrium observations curated in the same metabolomics study ( Bennett et al., 2009 ).
Proteomics data were taken from three sources to estimate cellular concentrations ( Schmidt et al., 2016;Taniguchi, 2010;Wi śniewski and Rakus, 2014 ). Again, only the dominant isozyme was included. All enzymes had at least two reported values.
Fluxomics data were not utilized as training data in this study. However, based on the work by Toya et al. (2010) , an output pyruvate production rate of 0.14 mMs −1 was selected as a constraint. Further, to capture the small albeit nontrivial effects of metabolite dilution due to cellular growth, we assume a cellular doubling time of one hour ( Woldringh et al., 1977 ). Finally, molecular weights are taken from EcoCyc ( Keseler et al., 2017 ). These molecular weights are not used directly in the dynamic model, but are utilized in scaling the optimization objective function.

A3. Determinacy analysis
The determinacy of a given parameter can be ascertained by inspecting F , the matrix that relates the basic parameters x to the observed parameter values b via b = F x . We define an undetermined parameter as one that is wholly uninformed by the training data; these are reflected by empty columns in F . A partially determined parameter is one that can only be explicitly defined in the context of fixing one or more other partially determined parameters. In the linearized form b = F x, this is equivalent to asking whether there exists a vector x such that F x = F (x + x ) , or equivalently F x = 0 , meaning that the perturbed parameter value vector x + x has the same fit as the original set of parameters x . Such vectors x are said to be in the null space of F , and therefore are composed by linear combinations of basis vectors for the null space.
This idea can be expressed compactly by forming a matrix N where the columns form a complete basis for the null space. Such a matrix can be found by singular value decomposition (SVD), or else by forming the null space orthogonal projection matrix where the symbol † indicates the Moore-Penrose pseudoinverse; in this case, N is an n × n matrix where n is the number of parameters (the size of x ). Regardless of how N is formed, we can use this matrix to find vectors η such that x = Nη for any x , where 0 = F x . If a row in N has a nonzero element, it means that the corresponding basic parameter can be perturbed in a fashion that does not change the overall fit to the data. Thus we call these parameters partially determined. The remaining parameters must then be fully determined.

A4. Parsimonious perturbation approach
The unit parsimonious perturbation vectors, used for exploring the parameter space in our optimization problem, can be calculated using the formula u = I − N ( AN ) † A d T † (A.12) where d is some desired perturbation direction (i.e. the basic or derived variable z = d T x that we wish to perturb, such that d T u = 1 ), A is a matrix which defines the basic or derived variables y = Ax that we wish to hold constant, N is a null space basis on the vector d T (see Eq. (A.11) for one way to obtain such a basis), and I is an identity matrix. Again, the symbol † indicates the Moore-Penrose pseudoinverse. The complete derivation of this expression can be found in the Supplementary Materials.
For the variables y (those which we generally wish to hold constant under a perturbation) we choose the energy differences associated with the forward and reverse maximum rates of reaction ( G max,f and G max,r ) as well as the binding energy differences for all substrates G b . We also include the energies from concentration G c for all metabolites with dynamic concentrations. Together these derived variables can fully parameterize the ODE system given in Eq. (1) . For the variables z = d T x (those which we wish to perturb) we choose all variables in y . We also include all basic variables such that the set of unit parsimonious perturbation vectors is full rank (discarding any redundant u vectors found in the process). These unit vectors are calculated once, and then reused throughout the optimization procedure.

A5. Objective formulation
Our parameter estimation objective, similar to that used by Jahan et al. (2016) , is thoroughly motivated and developed in the Supplementary Methods. In brief, the original problem may be stated as min x f (x ) (A.13) s.t. g(x ) = 0 (A.14) where f ( x ) is (again) our misfit cost term, and g ( x ) is a constraint term that is defined to be zero when we have obtained some nontrivial steady state. Since g ( x ) is nonlinear, g(x ) = 0 is difficult if not impossible to enforce, and thus we choose instead to solve the problem min x f (x ) + w · g(x ) (A.15) where we increase the weight w until we are satisfied that g ( x ) is sufficiently small.
We express f ( x ) as the sum of the absolute values of the differences in magnitude between parameter values and corresponding observations. The difference-in-magnitude is expressed as the difference between the logarithms of the estimated and observed values (or, equivalently, the logarithm of the ratio). Under our energetic reformulation, this simply becomes The difference-in-magnitude penalty is preferable because it is dimensionless (in contrast with, say, an absolute difference). We use the sum of absolute values, rather than the sum of squares, because the former tends to extract a few poor fits to the training data (i.e. points of frustration) rather than distributing the error more or less evenly across all observations.