BACKGROUND
high-throughput assays such as mass spectrometry are improving rapidly, which creates an opportunity for large scale in vivo measurements of the metabolome.
those in vivo data could enable estimation of kinetic parameters of metabolic reactions which are hard to estimate using in vitro data.
metabolic reactions are normally enzyme-catalyzed reactions, and quantitative estimates of their kinetic parameters could be very useful.
knowledge of kinetic parameters allows estimation of reaction rates directly from concentration measurements.
comparing the estimated kinetic parameters of a reaction in the wild type and mutant cells permits quantification of the effects of genetic variants, which may change the abundance or activity of a metabolic enzyme.
similarly, the effect of a drug that targets a particular enzyme could be estimated.
if parameters can be estimated for many reactions in a pathway, it would enable inference of the target of a mutation or drug treatment – if the estimates show that one enzyme is particularly strongly affected, that enzyme is probably the target.
finally, estimated parameters also allow estimation of maximum reaction rates, which can then be used as constraints to improve flux balance analysis  <cit> .
we explore the central problem of how to estimate the kinetic parameters of individual reactions using in vivo high-throughput measurements of metabolite concentrations and reaction rates at steady state, obtained by mass spectrometry or by nuclear magnetic resonance.
the method requires metabolite concentration and reaction rate data in multiple experiments under varying conditions.
for example, data could consist of several experiments obtained by perturbing the system through changes in nutrient media, drug treatment, or genetic alterations.
from such data, the kinetics of many individual reactions can potentially be estimated.
enzyme kinetic parameters have been measured for at least a century  <cit> .
the basic method involves mixing a measured amount of substrate and enzyme, and measuring the concentration of product at various points in time, creating a progress curve  <cit> .
in this setting, the experimenter has control over the initial concentrations of enzyme and substrate and thus can obtain relatively accurate measurements for concentrations.
although the experimental conditions are not at steady state, the mathematical formula for the kinetics can be simplified to the familiar michaelis menten kinetics by assuming that some elementary reactions are in near-equilibrium .
in contrast with an in vitro experiment, one major challenge with in vivo measurements of concentrations and reaction rates is the presence of significant error.
except for very low abundance metabolites, the errors are normally relative, meaning that they are proportional to the metabolite concentrations, instead of additive.
to quantify measurement precision with relative errors, experimentalists often use the coefficient of variation , which is calculated by dividing the standard deviation of peak area/height by the mean peak area/height .
methods such as least squares, which assume additive errors, are often not going to produce accurate estimates of parameters with relative errors.
because of such significant relative errors, it might not be reasonable to assume that errors are only in reaction rates as most of the in vitro enzyme kinetics methods assume.
relative errors in both concentrations and reaction rates need to be considered.
furthermore, many in vivo experiments are not time courses, so the data are assumed to be at steady-state.
another challenge with in vivo measurements is the difficulty of measuring enzyme and intermediate enzyme complex concentrations  <cit> , so these are typically unknown.
finally, control over metabolite concentrations in the cell is limited, so the range of experimental data points may be suboptimally distributed for accurate estimation of all parameters, making it difficult to estimate some parameters of a reaction.
a new estimation method, invest, standing for in vivo estimation, is described for estimating reaction parameters that addresses the specific challenges of in vivo data.
invest uses maximum likelihood estimation, based on a model where all measurements have relative errors.
as described, invest uses a family of reversible reaction mechanisms with multiple reactants and products with a single displacement mechanism.
it is not always possible to obtain data from the entire range of metabolite concentrations and reaction rates, so some parameters may not be identifiable.
invest estimates the standard deviations of parameter estimates using bootstrapping , so that the user can understand the range of errors for the estimates.
many methods for estimating kinetic parameters have been proposed, ranging from informal graphical plotting to sophisticated statistical non-linear regression methods.
however, none have addressed all of the problems of in vivo estimation discussed above.
many methods are based on the michaelis menten equation which are normally applied to irreversible single substrate, single product reactions.
standard graphical plotting methods, such as the double reciprocal plot  <cit>  and direct linear plot  <cit> , are not based on statistical estimation and yield unnecessarily inaccurate parameter estimates.
some more statistically-based methods deal with relative error or errors in all measurements – but not both.
specifically, weighted least squares  <cit>  is a general method often used in non-linear regression that can be applied to various kinds of reactions, however, it assumes the errors are additive and that only reaction rates have errors.
total least squares  <cit>  improves ordinary least squares by dealing with errors in all measurements, but the errors are still assumed to be additive.
raaijmakers’ maximum likelihood estimation method  <cit>  can deal with relative errors, but assumes that errors are in reaction rates only.
liebermeister et al.
<cit>  have developed a method that integrates knowledge from many sources, along with in vivo measurements, to estimate kinetic parameters using bayesian methods.
however, this method still assumes only additive errors and requires a lot of prior information about the parameters.
only invest deals with relative errors in all measurements as well as reversible reactions with multiple substrates and products.
a summary of existing methods appears in table  <dig>
in this paper, our goal is to focus on the specific problem of estimating kinetic parameters as accurately as possible, given realistic assumptions about data errors.
we discuss the formulation of invest, and evaluate the method on simulated data.
we show that invest works well on data with relative errors in all measurements.
we also demonstrate the application of invest and discuss the parameter identifiability issue.
methods
like most methods of kinetic parameter estimation, we assume that temperature and pressure are constant, so rate constants in mass action kinetic equations are constant, and the gibbs free energy of formation is constant.
we also assume that the measured system is at steady state, meaning that the time derivatives of metabolite concentrations and reaction rates are zero.
also, we assume that there are measurements of stable reactants and products of enzyme reactions, but not substrate-enzyme complexes, product-enzyme complexes and free enzyme concentrations, as they are generally difficult to measure experimentally.
it is assumed that metabolite concentrations are obtained by high-throughput methods, such as chromatography, mass spectroscopy, or nuclear magnetic resonance spectroscopy  <cit> .
for example, reasonably accurate concentration data can be obtained by mass spectroscopy with internal standards.
normally, average value of coefficient of variation for mass spectrometry below  <dig>  is considered as good measurements , and thus it is not unreasonable to expect such data to have a constant coefficient of variation  of  <dig> %.
we also assume that it is possible to obtain measurements of reaction rates.
for steady state reaction rate measurement, one widely used method is c <dig> labeling, which uses a cell culture at steady state in a medium with labeled-carbon substrates.
reaction rates can be determined by analyzing the labeling pattern of targeted metabolites from mass spectrometry  <cit> .
in addition, we assume that the gibb’s free energies of formation of metabolites are known, since these are used to compute the equilibrium constants  for enzymatic reactions.
single substrate and product reversible reactions
we use a standard simple but general reaction mechanism to represent most metabolic reversible reactions  <cit> .
this subsection considers single reactant/product case.
the more general case consisting of multiple reactants and multiple products will be discussed later.
the reaction is a three step process, namely binding, conversion and release:
a+e⇌k1k-1ae⇌k2k-2be⇌k3k-3b+e
where a is the reactant, b is the product, e is the free enzyme, ae and be are the intermediate complexes, and ki and k−i are reaction rate constants for i∈{ <dig> ,3}.
assuming the reaction is at steady state, an equation for the reaction rate can be written as:
v=keq−c1+c2+c <dig>
where keq=k1k2k3k−1k−2k− <dig> is an equilibrium constant, obtained from the standard gibbs free energy of formation of the reactants and products.
c <dig> is
k2k3k−1k−2k−3+k3k−2k−3+1k−3/etot,
c <dig> is
k1k2k−1k−2k−3+k1k3k−1k−2k−3+k1k−1k−3/,
c <dig> is
1k−2+1k−1+k2k−1k−2/,  and , the total enzyme, is ++.
if keq is very large and the reversible reactions’ rate constants  are small, c <dig> can be neglected and the rate eq.
<dig> can be reduced to standard irreversible michaelis menten equation.
this rate equation can be derived from the ordinary differential equations for mass action kinetics of a reaction , by setting the derivatives of the concentrations of all chemical species to zero  and solving for .
the detailed derivation and calculation for the steady state equation and equilibrium constant are presented in additional files  <dig> and  <dig>
parameter estimation by maximum likelihood for single substrate/product reversible reaction
the invest method estimates the parameters of kinetic rate eq.
using maximum likelihood, assuming relative error in all measurements.
parameters are estimated from a set of n experiments, each with data values for ai , bi , vi , for experiment i.
each data value has some known relative error.
specifically, we have ai=ai0εa, bi=bi0εb and vi=vi0εv, where ai <dig>  bi <dig>  and vi <dig> are latent variables representing the data values without measurement error, multiplied by a normally distributed error with mean  <dig> and standard deviation σ: εx∼n <dig> σx <dig> .
the likelihood function is:
l=f
since each data acquisition can be carried out independently  <cit> , errors in a, b and v can be assumed to be independent of c <dig> c <dig> and c <dig> and each other, the likelihood function can be written as
f=∏f∏f∏f
the distribution of ai is
nai <dig> ai02σa2=12πai02σa2exp−22ai02σa <dig>  the distributions of the other data values are similar.
the parameters that maximize the likelihood also maximize the log of the likelihood, which is
log=∑i=1n−logai0σa2π+∑i=1n−22ai02σa2+∑i=1n−logbi0σb2π+∑i=1n−22bi02σb2+∑i=1n−logvi0σv2π+∑i=1n−22vi02σv <dig>
negating the log likelihood and dropping constant factors yields an objective function to minimize, subject to the constraints of eq.
<dig>
min∑i=1n+log+log)+12σa2∑i=1naiai0−12+12σb2∑i=1nbibi0−12+12σv2∑i=1nvivi0−12s.t.vi0=keqai0−bi0c1+c2ai0+c3bi <dig> wherei= <dig> ,⋯,n  where all the ai, bi and vi are experimental measurements, all the relative errors σ are known and ai <dig>  bi <dig>  vi <dig> are latent variables, and c <dig>  c <dig> and c <dig> are the parameters to be estimated by solving the optimization problem.
in the implementation, this is simplified to an unconstrained optimization problem by substituting the right-hand side of eq.
<dig> for vi <dig>
generalization to multiple substrates and products
for reactions with multiple substrates and products, there are two possible mechanisms, namely single-displacement and double-displacement.
for single-displacement reactions, the order of substrates binding to the enzyme can be random or ordered.
those two type of reactions can be approximated by following reaction  <cit> :
a1+a2+⋯+am+e⇌k1k-1a1a2⋯ame⇌k2k-2b1b2⋯bpe⇌k3k-3b1+b2+⋯+bp+e  where m is the number of reactants and p is the number of products in this reaction.
a steady state equation can be derived as in the single reactant/product case:
v=keq∏j=1m−∏j=1pc1+c2∏j=1m+c3∏j=1p
where c <dig>  c <dig>  c <dig>  keq, and etot are as before.
the derivation of the objective function to minimize in order to find the parameters that maximize the likelihood is a straightforward generalization of the single substrate/product case.
min∑i=1n∑j=1mlog+∑i=1n∑j=1plog+∑i=1nlog+12σa2∑i=1n∑j=1maijaij0−12+12σb2∑i=1n∑j=1pbijbij0−12+12σv2∑i=1nvivi0− <dig>  which is maximized subject to the constraints of eq.
<dig>
in the implementation, this can also be simplified to an unconstrained optimization problem by substituting the right-hand side of eq.
<dig> for vi <dig>
parameter identifiability
it is sometimes not possible to obtain in vivo data whose values are well enough distributed to estimate all parameters accurately.
in this section, we characterize some cases when parameters cannot be accurately estimated.
from eq.
, it is clear that when one term in the denominator is much smaller than the others, v is relatively insensitive to the corresponding parameter.
for example, if c <dig> c2a≫c3b, then eq.
<dig> will be approximately
v=keq−c1+c <dig>
so changes in c <dig> will have little effect on v. more importantly, changes in data values resulting from erroneous estimates of c <dig> will be small relative to the noise in the data, so estimates of c <dig> tend to have large errors.
similarly, estimates of c <dig> tend to have large errors when c2a+c3b≫c <dig> and estimates of c <dig> have large errors when c1+c3b≫c2a.
for illustration, consider the simpler case when keq is very large and the rate eq.
can be approximated by the standard michaelis menten equation.
in fig.
<dig>  two data sets derived from the same actual parameters have large ai, so the vi values lie near the maximum value of the curve.
we call this region as saturation region since reaction rates asymptotically approach a maximum level, and additional increases in the substrate concentration do not lead to an increase in the reaction rates.
in this case, c <dig>  which determines the maximum value, is the only parameter that affects the curve fit, so estimates of c <dig> from both data sets have large errors.
in fig.
<dig>  all of the substrate concentration ai values are small, so the points lie near the region where the curve is increasing linearly.
we call this region as linear region since reaction rates increase in almost a linear fashion with increasing substrate concentrations.
the slope in this region is determined by c <dig> almost independently of c <dig> so estimates of c <dig> have large errors.
fig.
<dig> identifiability issue in two parameter case.
when data points are not well-distributed, parameter identification can be difficult.
this shows the curve for parameters predicted from two possible data sets, one with points gathered in the saturation region  in  and in the other in the linear region  in
estimates of the accuracy of parameter estimates must be obtained using the available data.
invest uses bootstrapping to estimate the variance of the parameter estimates.
bootstrap estimation of standard error
the c parameter estimates can vary widely in accuracy, depending on the experimental data.
bootstrapping  <cit>  is used to estimate the relative standard errors and bias of the parameter estimates, so users can tell whether the parameter estimation is good or not.
let ĉ be the estimate from the data, and ci^∗ be the estimate from a bootstrap sample.
a typical recommendation is to use n=n <dig> bootstrap samples for n experimental measurements  <cit> .
the bootstrap estimation of standard errors is calculated from seb=1n∑ <dig> and bias estimation is calculated by bias=1n∑ci^∗−ĉ <cit> .
as the c parameters have a large range of possible values, it is more appropriate to use relative errors and relative bias to describe the estimate.
the relative standard error is calculated by seb/ĉ and the relative bias is calculated by bias/ĉ.
estimation of total enzyme change
estimating kinetic parameters can be useful for identifying the effects of genetic changes or drug treatments that target metabolic enzymes.
the total concentration of the enzyme in the cell may change because of changes in gene expression or loss of function in one or more copies of the gene coding for the enzyme, or the activity may change because of changes in the protein sequence or post translational modifications.
estimating these changes for specific enzymes in each sample can help identify the target of a mutation or drug , and may be useful for estimating the impact of such a change on flux through a network.
since each of the kinetic parameters ci is of the form ci′/etot, where ci′ is independent of the enzyme concentration, etot can be estimated from the ratio
etotwtetotmt=cimtciwt  where ciwt and cimt are corresponding ci parameters  for wild type and mutant  samples.
note that it is possible to obtain a reliable estimate for etot whenever there are reliable estimates for one of the three parameters in both samples.
RESULTS
evaluate invest using simulated data
we evaluate the parameter estimation method on simulated data.
for matlab code for reproducing the results of this work, please refer to  <cit> .
the simulations were carried out in matlab on a laptop computer with an intel core i5-4200u  <dig>  ghz processor and  <dig> gb installed memory.
many reactions in metabolic pathways have multiple substrates and products and are reversible reactions.
the simulation is based on the reaction acetylornithine aminotransferase from saccharomyces cerevisiae arginine biosynthesis pathway with arg <dig>  <cit> .
kinetic parameters and the total enzyme concentration are not available, and thus we use some heuristic numbers for them.
the experimental data are chosen to be well-distributed, since poorly distributed data would guarantee inaccurate parameter estimates even for the best possible estimation method.
the reaction is:
acglu-sa+glu⇌acorn+2−oxoglutarate  abbreviations  <cit> : acglu-sa, n-acetyl-glutamate-semialdehyde; glu, l-glutamate; acorn, n-acetyl- ornithine.
the standard gibbs free energy of formation for the metabolites are taken from metacyc database  <cit> , and are provided in the additional file  <dig>  the standard gibbs free energy of formation can be used to compute keq= <dig> , and, assuming etot= <dig> m, the c parameters are c1= <dig> , c2= <dig>  and c3= <dig> .
to characterize the amount of data for effective use of invest, we evaluated the accuracy of parameter estimates for varying numbers of simulated experiments.
data sets of  <dig>   <dig> and  <dig> experiments were generated by choosing values for substrate and product concentrations and computing v exactly for each choice based on eq.
<dig>  relative errors were introduced by multiplying a random value from the normal distribution of n. a value of  <dig>  was used for σ for metabolites, and σv of  <dig>  was used for reaction rates.
for each number of experiments,  <dig>  simulated data sets were generated, the c parameters were estimated using invest, and the mean and standard error were calculated.
the results are shown in table  <dig>  with increasing sample size, the relative standard errors and bias in the estimates are improved.
it is evident that the results for sample size of  <dig> and  <dig> are quite accurate with relative standard error near  <dig> % and very small relative bias.
twenty to thirty samples seems to be a reasonable sample size to choose for accurate estimations.v= <dig> .
results are based on  <dig>  simulated data sets.
“n” is the number of experiments.
“avg est” is the average value of the estimates.
“rel se” is the relative standard error, and “rel bias” is the relative bias
c
c
c
c
c
c
c
c
c
second, we consider the effect of greater error in reaction rate estimates, with σv= <dig> .
the results are shown in table  <dig>  the relative standard errors increase, but are still below  <dig> %.
the relative bias values are also low.
this shows that invest is robust to different measurement errors.c parameter estimates for acetylornithine aminotransferase when σ
v= <dig> .
results are based on  <dig> simulated data sets of  <dig> experiments, each
c
1
c
2
c
3
it is also possible to evaluate the accuracy of estimates when there is only one data set  available, as would be the case in normal use of invest in practice.
the bootstrap method is used to estimate relative standard errors in parameter estimates.
to evaluate the bootstrap method, we generated a single data set of  <dig> experiments as the input data for parameter estimation and randomly subsampled the  <dig> data points  <dig>  times.
each bootstrap subsample simulation took around  <dig> sec.
the estimates for σv= <dig>  and σv= <dig>  are shown in tables  <dig> and  <dig> respectively.
as expected, the bootstrap estimates are very similar to the previous estimates from  <dig>  simulated data sets.c parameter estimates for acetylornithine aminotransferase when σ
v= <dig> .
estimates are from a single simulated data set of  <dig> experiments.
the bootstrap method was used to estimate relative standard error  and relative bias
c
1
c
2
c
3
v= <dig> .
estimates are from a single simulated data set of  <dig> experiments.
the bootstrap method was used to estimate relative standard error  and relative bias
c
1
c
2
c
3
comparison of invest with prior methods
most current methods produce optimal estimates only when errors are additive and when errors occur only in reaction rate measurements.
these assumptions are generally not true with in vivo data.
in this subsection, we compare invest to some existing methods and show that invest produces better estimates when data have relative errors in all measurements.
as some of the existing methods only work on irreversible enzymatic reactions, we use the two parameter case of eq.
<dig> for comparison.
in this case, there are two parameters to be estimated, namely c <dig> and c <dig>
v=keqc1+c <dig>
we first simulate the data with relative errors to both substrate a and reaction rate v, and second apply invest and prior methods to obtain estimates for the michaelis menten like curve.
one thousand simulated data sets of  <dig> experiments each are used.
the results are summarized in table  <dig>  invest has superior performance in the estimates and relative standard errors.1= <dig> , true c
2= <dig> .
data have relative errors in all variables.
results are based on  <dig>  simulated data sets of  <dig> experiments, each.
“avg est” is the average value of the estimates.
“rel se” is the relative standard error
1
2
1
2
predicting total enzyme concentration change
as noted above, the relative difference in etot between wild type and mutant or drug-treated samples can be estimated from the estimate of any of the ci parameters from two sets of experiments.
etotwtetotmt=cimtciwt.
we illustrate estimation of etot change using the arg <dig> reaction.
for the wild type samples, the total enzyme concentration is etotwt=1m, and for the mutant/drug treated samples, the total enzyme concentration is etotmt= <dig> m. results of the wild-type estimate appear in the previous section.
additional data for the mutant were generated as above based on the c parameter values of mutant/drug treated sample and  <dig>  simulated data sets are used.
the estimates for mutant/drug treated sample are shown in table  <dig> c parameter estimates for acetylornithine aminotransferase from mutant/drug treated sample.
results are based on  <dig>  simulated data sets
c
1
c
2
c
3
to obtain the prediction of total enzyme change, we take cimt/ciwt.
the results are shown in table  <dig> e
tot change prediction based on  <dig>  simulated data sets
etotwtetotmt=c1mtc1wt
etotwtetotmt=c2mtc2wt
etotwtetotmt=c3mtc3wt
since any of the ci parameters can be used to estimate the change in etot, the one that gives minimum standard error, c <dig>  was chosen.
this also demonstrates that even though sometimes identifiability issues can occur and some parameters cannot be estimated, our method could still be very useful if one parameter can be estimated accurately.
discussion
this work is intended to be a first step towards estimating parameters for reactions in large metabolic networks in vivo.
in vivo estimation will need to be based on data that have relatively large relative errors in all measured parameters, and will have to deal with a variety of reaction kinetics, including reactions that are reversible and have multiple substrates and/or products.
although measurement and estimation of enzyme kinetics has been studied for many decades, there is no single existing estimation method that addresses all of these issues.
we have proposed a maximum likelihood approach to estimate kinetic parameters using nonlinear optimization, with estimates on the standard error and bias of the results using the bootstrap.
simulations show that invest produces accurate estimates for realistic high-throughput metabolomics data.
for example, with 20– <dig> samples with coefficients of  <dig> % in metabolite concentrations and  <dig> % in reaction rate estimates, estimates have a relative standard error of less than  <dig> %.
collecting data of this quality would be technically difficult, but is within the current state of the art.
an advantage of the method is that it estimates each set of reaction parameters independently.
if measurements are not available for some metabolites, it can still estimate parameters for those reactions for which the data include all substrates and products.
solving the problem of in vivo parameter estimation in its full generality will require meeting a number of additional challenges.
some reactions have more complex kinetics than those we consider, especially various kinds of inhibition.
when the inhibiting metabolite and mechanism of inhibition are known, the approach described here can probably be generalized to accommodate the inhibition mechanism in our future work.
otherwise, a process of model selection may be necessary, where competing models are estimated and the quality of the results compared, with appropriate adjustments for model complexity.
in addition, it will be necessary to deal with the kinetics of transport reactions, and to take account of different compartments in the cell.
parameter identifiability is a difficult issue in in vivo estimation.
we have shown that accurate estimates of all parameters require data that is well-distributed over the kinetics curve, but such data will not often be obtainable for several reasons.
experimental data must be obtained by perturbing metabolites and fluxes, for example, by adjusting nutrient media, testing mutants, and targeting reactions with drugs.
first, accurate estimation may require non-physiological concentrations of metabolites – estimating c <dig> for a reaction that is nearly irreversible being an example.
more generally, there is usually inadequate controllability of metabolite concentrations and reaction fluxes to obtain the experimental values needed for accurate estimation, for many reasons including concentrations are toxic or inadequate to sustain life, and rate-limiting reactions that make high fluxes in other reactions impossible to obtain.
since we can’t estimate everything accurately, it is important to produce estimates of the standard errors of parameter estimates, so we can tell which ones are meaningful.
also, as we note above, if some but not all parameters of a reaction can be estimated accurately, the results still may be useful.
for example, it is possible to estimate the total concentration or relative activity of an enzyme in wild-type vs. mutant cells when only one of the kinetic parameters is accurately estimated.
CONCLUSIONS
in conclusion, a new method, invest, is developed for estimating reaction kinetic parameters in metabolic networks that addresses the specific challenges of in vivo data.
invest uses maximum likelihood estimation based on models where all measurements have potentially relative errors.
it can be applied to multiple substrate/product reversible enzymatic reactions with a generalized single displacement mechanism.
because it is not always possible to obtain good data covering full range of possible metabolite concentrations and reaction rates, certain parameters may be non-identifiable.
invest uses bootstrap to estimate the standard errors of parameter estimations that can tell which estimates are reliable.
invest enables the estimation of reaction rates directly from concentration measurements.
also, comparing the estimated kinetic parameters of a reaction in the wild type and mutant cells can quantify enzyme abundance or activity change due to genetic variants.
the same method can also be used to measure the effect of a drug that targets a particular enzyme.
moreover, estimated parameters can be used to estimate maximum reaction rates, which could be used as constraints to improve flux-balance analysis.
additional files
additional file  <dig>
experimental data support on relative error model.
figure s <dig> noise errors in high-throughput metabolomic data tend to be relative.
the plot shows the empirical standard deviation vs. mean of metabolite concentrations in a publicly available mass spectrometry data set of  <dig> human urine samples  <cit> .
each sample has  <dig> technical replicates, which were used to calculate the standard deviation and mean of metabolite concentrations.
the data for “peak 105” were chosen because the chromatographic peak appears in all three replicates of the sample and the measurements cover a wide range of concentrations across different samples.
low concentrations are omitted because they are highly inaccurate due to background noise.
there is a linear relationship  between standard deviation and concentration mean, showing that errors are proportional to measured concentration.
additional file  <dig>
derivation for steady state rate equation.
this file provides a detailed derivation for steady state rate equation of a single reactant and single product reversible metabolic reaction presented in methods section.
additional file  <dig>
equilibrium constant
k
eq
.
the equilibrium constant k
eq is assumed to be a known constant.
this file provides the calculation of equilibrium constant k
eq based on standard gibbs free energies of formation.
additional file  <dig>
standard gibbs free energy of formation metacyc.
this file provides standard gibbs free energy of formation taken from metacyc database  <cit>  for metabolites used in the simulation example in results section.
availability of data and materials
matlab code for reproducing the results of this work is available at stanford digital repository, http://purl.stanford.edu/bg158sn <dig>
competing interests
the authors declare that they have no competing interests.
author’s contributions
wz and dld defined the problem.
wz did the mathematical derivations, implemented the method, and performed the simulations.
rk suggested the maximum likelihood approach and worked out key ideas about the approach.
dld oversaw the work.
the authors collaborated on the writing and all authors have read and approved the final manuscript.

 high-throughput assays such as mass spectrometry are improving rapidly, which creates an opportunity for large scale in vivo measurements of the metabolome. those in vivo data could enable estimation of kinetic parameters of metabolic reactions which are hard to estimate using in vitro data. metabolic reactions are normally enzyme-catalyzed reactions, and quantitative estimates of their kinetic parameters could be very useful. knowledge of kinetic parameters allows estimation of reaction rates directly from concentration measurements. comparing the estimated kinetic parameters of a reaction in the wild type and mutant cells permits quantification of the effects of genetic variants, which may change the abundance or activity of a metabolic enzyme. similarly, the effect of a drug that targets a particular enzyme could be estimated. if parameters can be estimated for many reactions in a pathway, it would enable inference of the target of a mutation or drug treatment – if the estimates show that one enzyme is particularly strongly affected, that enzyme is probably the target. finally, estimated parameters also allow estimation of maximum reaction rates, which can then be used as constraints to improve flux balance analysis. we explore the central problem of how to estimate the kinetic parameters of individual reactions using in vivo high-throughput measurements of metabolite concentrations and reaction rates at steady state, obtained by mass spectrometry or by nuclear magnetic resonance. the method requires metabolite concentration and reaction rate data in multiple experiments under varying conditions. for example, data could consist of several experiments obtained by perturbing the system through changes in nutrient media, drug treatment, or genetic alterations. from such data, the kinetics of many individual reactions can potentially be estimated. enzyme kinetic parameters have been measured for at least a century. the basic method involves mixing a measured amount of substrate and enzyme, and measuring the concentration of product at various points in time, creating a progress curve. in this setting, the experimenter has control over the initial concentrations of enzyme and substrate and thus can obtain relatively accurate measurements for concentrations. although the experimental conditions are not at steady state, the mathematical formula for the kinetics can be simplified to the familiar michaelis menten kinetics by assuming that some elementary reactions are in near-equilibrium. in contrast with an in vitro experiment, one major challenge with in vivo measurements of concentrations and reaction rates is the presence of significant error. except for very low abundance metabolites, the errors are normally relative, meaning that they are proportional to the metabolite concentrations, instead of additive. to quantify measurement precision with relative errors, experimentalists often use the coefficient of variation, which is calculated by dividing the standard deviation of peak area/height by the mean peak area/height. methods such as least squares, which assume additive errors, are often not going to produce accurate estimates of parameters with relative errors. because of such significant relative errors, it might not be reasonable to assume that errors are only in reaction rates as most of the in vitro enzyme kinetics methods assume. relative errors in both concentrations and reaction rates need to be considered. furthermore, many in vivo experiments are not time courses, so the data are assumed to be at steady-state. another challenge with in vivo measurements is the difficulty of measuring enzyme and intermediate enzyme complex concentrations, so these are typically unknown. finally, control over metabolite concentrations in the cell is limited, so the range of experimental data points may be suboptimally distributed for accurate estimation of all parameters, making it difficult to estimate some parameters of a reaction. a new estimation method, invest, standing for in vivo estimation, is described for estimating reaction parameters that addresses the specific challenges of in vivo data. invest uses maximum likelihood estimation, based on a model where all measurements have relative errors. as described, invest uses a family of reversible reaction mechanisms with multiple reactants and products with a single displacement mechanism. it is not always possible to obtain data from the entire range of metabolite concentrations and reaction rates, so some parameters may not be identifiable. invest estimates the standard deviations of parameter estimates using bootstrapping, so that the user can understand the range of errors for the estimates. many methods for estimating kinetic parameters have been proposed, ranging from informal graphical plotting to sophisticated statistical non-linear regression methods. however, none have addressed all of the problems of in vivo estimation discussed above. many methods are based on the michaelis menten equation which are normally applied to irreversible single substrate, single product reactions. standard graphical plotting methods, such as the double reciprocal plot and direct linear plot, are not based on statistical estimation and yield unnecessarily inaccurate parameter estimates. some more statistically-based methods deal with relative error or errors in all measurements – but not both. specifically, weighted least squares is a general method often used in non-linear regression that can be applied to various kinds of reactions, however, it assumes the errors are additive and that only reaction rates have errors. total least squares improves ordinary least squares by dealing with errors in all measurements, but the errors are still assumed to be additive. raaijmakers’ maximum likelihood estimation method can deal with relative errors, but assumes that errors are in reaction rates only. liebermeister et al. have developed a method that integrates knowledge from many sources, along with in vivo measurements, to estimate kinetic parameters using bayesian methods. however, this method still assumes only additive errors and requires a lot of prior information about the parameters. only invest deals with relative errors in all measurements as well as reversible reactions with multiple substrates and products. a summary of existing methods appears in table in this paper, our goal is to focus on the specific problem of estimating kinetic parameters as accurately as possible, given realistic assumptions about data errors. we discuss the formulation of invest, and evaluate the method on simulated data. we show that invest works well on data with relative errors in all measurements. we also demonstrate the application of invest and discuss the parameter identifiability issue. methods like most methods of kinetic parameter estimation, we assume that temperature and pressure are constant, so rate constants in mass action kinetic equations are constant, and the gibbs free energy of formation is constant. we also assume that the measured system is at steady state, meaning that the time derivatives of metabolite concentrations and reaction rates are zero. also, we assume that there are measurements of stable reactants and products of enzyme reactions, but not substrate-enzyme complexes, product-enzyme complexes and free enzyme concentrations, as they are generally difficult to measure experimentally. it is assumed that metabolite concentrations are obtained by high-throughput methods, such as chromatography, mass spectroscopy, or nuclear magnetic resonance spectroscopy. for example, reasonably accurate concentration data can be obtained by mass spectroscopy with internal standards. normally, average value of coefficient of variation for mass spectrometry below is considered as good measurements, and thus it is not unreasonable to expect such data to have a constant coefficient of variation of %. we also assume that it is possible to obtain measurements of reaction rates. for steady state reaction rate measurement, one widely used method is c labeling, which uses a cell culture at steady state in a medium with labeled-carbon substrates. reaction rates can be determined by analyzing the labeling pattern of targeted metabolites from mass spectrometry. in addition, we assume that the gibb’s free energies of formation of metabolites are known, since these are used to compute the equilibrium constants for enzymatic reactions. single substrate and product reversible reactions we use a standard simple but general reaction mechanism to represent most metabolic reversible reactions. this subsection considers single reactant/product case. the more general case consisting of multiple reactants and multiple products will be discussed later. the reaction is a three step process, namely binding, conversion and release: a+e⇌k1k-1ae⇌k2k-2be⇌k3k-3b+e where a is the reactant, b is the product, e is the free enzyme, ae and be are the intermediate complexes, and ki and k−i are reaction rate constants for i∈{,3}. assuming the reaction is at steady state, an equation for the reaction rate can be written as: v=keq−c1+c2+c where keq=k1k2k3k−1k−2k− is an equilibrium constant, obtained from the standard gibbs free energy of formation of the reactants and products. c is k2k3k−1k−2k−3+k3k−2k−3+1k−3/etot, c is k1k2k−1k−2k−3+k1k3k−1k−2k−3+k1k−1k−3/, c is 1k−2+1k−1+k2k−1k−2/, and, the total enzyme, is ++. if keq is very large and the reversible reactions’ rate constants are small, c can be neglected and the rate eq. can be reduced to standard irreversible michaelis menten equation. this rate equation can be derived from the ordinary differential equations for mass action kinetics of a reaction, by setting the derivatives of the concentrations of all chemical species to zero and solving for. the detailed derivation and calculation for the steady state equation and equilibrium constant are presented in additional files and parameter estimation by maximum likelihood for single substrate/product reversible reaction the invest method estimates the parameters of kinetic rate eq. using maximum likelihood, assuming relative error in all measurements. parameters are estimated from a set of n experiments, each with data values for ai, bi, vi, for experiment i. each data value has some known relative error. specifically, we have ai=ai0εa, bi=bi0εb and vi=vi0εv, where ai bi and vi are latent variables representing the data values without measurement error, multiplied by a normally distributed error with mean and standard deviation σ: εx∼n σx. the likelihood function is: l=f since each data acquisition can be carried out independently, errors in a, b and v can be assumed to be independent of c c and c and each other, the likelihood function can be written as f=∏f∏f∏f the distribution of ai is nai ai02σa2=12πai02σa2exp−22ai02σa the distributions of the other data values are similar. the parameters that maximize the likelihood also maximize the log of the likelihood, which is log=∑i=1n−logai0σa2π+∑i=1n−22ai02σa2+∑i=1n−logbi0σb2π+∑i=1n−22bi02σb2+∑i=1n−logvi0σv2π+∑i=1n−22vi02σv negating the log likelihood and dropping constant factors yields an objective function to minimize, subject to the constraints of eq. min∑i=1n+log+log)+12σa2∑i=1naiai0−12+12σb2∑i=1nbibi0−12+12σv2∑i=1nvivi0−12s.t.vi0=keqai0−bi0c1+c2ai0+c3bi wherei=,⋯,n where all the ai, bi and vi are experimental measurements, all the relative errors σ are known and ai bi vi are latent variables, and c c and c are the parameters to be estimated by solving the optimization problem. in the implementation, this is simplified to an unconstrained optimization problem by substituting the right-hand side of eq. for vi generalization to multiple substrates and products for reactions with multiple substrates and products, there are two possible mechanisms, namely single-displacement and double-displacement. for single-displacement reactions, the order of substrates binding to the enzyme can be random or ordered. those two type of reactions can be approximated by following reaction: a1+a2+⋯+am+e⇌k1k-1a1a2⋯ame⇌k2k-2b1b2⋯bpe⇌k3k-3b1+b2+⋯+bp+e where m is the number of reactants and p is the number of products in this reaction. a steady state equation can be derived as in the single reactant/product case: v=keq∏j=1m−∏j=1pc1+c2∏j=1m+c3∏j=1p where c c c keq, and etot are as before. the derivation of the objective function to minimize in order to find the parameters that maximize the likelihood is a straightforward generalization of the single substrate/product case. min∑i=1n∑j=1mlog+∑i=1n∑j=1plog+∑i=1nlog+12σa2∑i=1n∑j=1maijaij0−12+12σb2∑i=1n∑j=1pbijbij0−12+12σv2∑i=1nvivi0− which is maximized subject to the constraints of eq. in the implementation, this can also be simplified to an unconstrained optimization problem by substituting the right-hand side of eq. for vi parameter identifiability it is sometimes not possible to obtain in vivo data whose values are well enough distributed to estimate all parameters accurately. in this section, we characterize some cases when parameters cannot be accurately estimated. from eq., it is clear that when one term in the denominator is much smaller than the others, v is relatively insensitive to the corresponding parameter. for example, if c c2a≫c3b, then eq. will be approximately v=keq−c1+c so changes in c will have little effect on v. more importantly, changes in data values resulting from erroneous estimates of c will be small relative to the noise in the data, so estimates of c tend to have large errors. similarly, estimates of c tend to have large errors when c2a+c3b≫c and estimates of c have large errors when c1+c3b≫c2a. for illustration, consider the simpler case when keq is very large and the rate eq. can be approximated by the standard michaelis menten equation. in fig. two data sets derived from the same actual parameters have large ai, so the vi values lie near the maximum value of the curve. we call this region as saturation region since reaction rates asymptotically approach a maximum level, and additional increases in the substrate concentration do not lead to an increase in the reaction rates. in this case, c which determines the maximum value, is the only parameter that affects the curve fit, so estimates of c from both data sets have large errors. in fig. all of the substrate concentration ai values are small, so the points lie near the region where the curve is increasing linearly. we call this region as linear region since reaction rates increase in almost a linear fashion with increasing substrate concentrations. the slope in this region is determined by c almost independently of c so estimates of c have large errors. fig. identifiability issue in two parameter case. when data points are not well-distributed, parameter identification can be difficult. this shows the curve for parameters predicted from two possible data sets, one with points gathered in the saturation region in and in the other in the linear region in estimates of the accuracy of parameter estimates must be obtained using the available data. invest uses bootstrapping to estimate the variance of the parameter estimates. bootstrap estimation of standard error the c parameter estimates can vary widely in accuracy, depending on the experimental data. bootstrapping is used to estimate the relative standard errors and bias of the parameter estimates, so users can tell whether the parameter estimation is good or not. let ĉ be the estimate from the data, and ci^∗ be the estimate from a bootstrap sample. a typical recommendation is to use n=n bootstrap samples for n experimental measurements. the bootstrap estimation of standard errors is calculated from seb=1n∑ and bias estimation is calculated by bias=1n∑ci^∗−ĉ. as the c parameters have a large range of possible values, it is more appropriate to use relative errors and relative bias to describe the estimate. the relative standard error is calculated by seb/ĉ and the relative bias is calculated by bias/ĉ. estimation of total enzyme change estimating kinetic parameters can be useful for identifying the effects of genetic changes or drug treatments that target metabolic enzymes. the total concentration of the enzyme in the cell may change because of changes in gene expression or loss of function in one or more copies of the gene coding for the enzyme, or the activity may change because of changes in the protein sequence or post translational modifications. estimating these changes for specific enzymes in each sample can help identify the target of a mutation or drug, and may be useful for estimating the impact of such a change on flux through a network. since each of the kinetic parameters ci is of the form ci′/etot, where ci′ is independent of the enzyme concentration, etot can be estimated from the ratio etotwtetotmt=cimtciwt where ciwt and cimt are corresponding ci parameters for wild type and mutant samples. note that it is possible to obtain a reliable estimate for etot whenever there are reliable estimates for one of the three parameters in both samples.  evaluate invest using simulated data we evaluate the parameter estimation method on simulated data. for matlab code for reproducing the results of this work, please refer to. the simulations were carried out in matlab on a laptop computer with an intel core i5-4200u ghz processor and gb installed memory. many reactions in metabolic pathways have multiple substrates and products and are reversible reactions. the simulation is based on the reaction acetylornithine aminotransferase from saccharomyces cerevisiae arginine biosynthesis pathway with arg. kinetic parameters and the total enzyme concentration are not available, and thus we use some heuristic numbers for them. the experimental data are chosen to be well-distributed, since poorly distributed data would guarantee inaccurate parameter estimates even for the best possible estimation method. the reaction is: acglu-sa+glu⇌acorn+2−oxoglutarate abbreviations: acglu-sa, n-acetyl-glutamate-semialdehyde; glu, l-glutamate; acorn, n-acetyl- ornithine. the standard gibbs free energy of formation for the metabolites are taken from metacyc database, and are provided in the additional file the standard gibbs free energy of formation can be used to compute keq=, and, assuming etot= m, the c parameters are c1=, c2= and c3=. to characterize the amount of data for effective use of invest, we evaluated the accuracy of parameter estimates for varying numbers of simulated experiments. data sets of and experiments were generated by choosing values for substrate and product concentrations and computing v exactly for each choice based on eq. relative errors were introduced by multiplying a random value from the normal distribution of n. a value of was used for σ for metabolites, and σv of was used for reaction rates. for each number of experiments, simulated data sets were generated, the c parameters were estimated using invest, and the mean and standard error were calculated. the results are shown in table with increasing sample size, the relative standard errors and bias in the estimates are improved. it is evident that the results for sample size of and are quite accurate with relative standard error near % and very small relative bias. twenty to thirty samples seems to be a reasonable sample size to choose for accurate estimations.v=. results are based on simulated data sets. “n” is the number of experiments. “avg est” is the average value of the estimates. “rel se” is the relative standard error, and “rel bias” is the relative bias c c c c c c c c c second, we consider the effect of greater error in reaction rate estimates, with σv=. the results are shown in table the relative standard errors increase, but are still below %. the relative bias values are also low. this shows that invest is robust to different measurement errors.c parameter estimates for acetylornithine aminotransferase when σ v=. results are based on simulated data sets of experiments, each c 1 c 2 c 3 it is also possible to evaluate the accuracy of estimates when there is only one data set available, as would be the case in normal use of invest in practice. the bootstrap method is used to estimate relative standard errors in parameter estimates. to evaluate the bootstrap method, we generated a single data set of experiments as the input data for parameter estimation and randomly subsampled the data points times. each bootstrap subsample simulation took around sec. the estimates for σv= and σv= are shown in tables and respectively. as expected, the bootstrap estimates are very similar to the previous estimates from simulated data sets.c parameter estimates for acetylornithine aminotransferase when σ v=. estimates are from a single simulated data set of experiments. the bootstrap method was used to estimate relative standard error and relative bias c 1 c 2 c 3 v=. estimates are from a single simulated data set of experiments. the bootstrap method was used to estimate relative standard error and relative bias c 1 c 2 c 3 comparison of invest with prior methods most current methods produce optimal estimates only when errors are additive and when errors occur only in reaction rate measurements. these assumptions are generally not true with in vivo data. in this subsection, we compare invest to some existing methods and show that invest produces better estimates when data have relative errors in all measurements. as some of the existing methods only work on irreversible enzymatic reactions, we use the two parameter case of eq. for comparison. in this case, there are two parameters to be estimated, namely c and c v=keqc1+c we first simulate the data with relative errors to both substrate a and reaction rate v, and second apply invest and prior methods to obtain estimates for the michaelis menten like curve. one thousand simulated data sets of experiments each are used. the results are summarized in table invest has superior performance in the estimates and relative standard errors.1=, true c 2=. data have relative errors in all variables. results are based on simulated data sets of experiments, each. “avg est” is the average value of the estimates. “rel se” is the relative standard error 1 2 1 2 predicting total enzyme concentration change as noted above, the relative difference in etot between wild type and mutant or drug-treated samples can be estimated from the estimate of any of the ci parameters from two sets of experiments. etotwtetotmt=cimtciwt. we illustrate estimation of etot change using the arg reaction. for the wild type samples, the total enzyme concentration is etotwt=1m, and for the mutant/drug treated samples, the total enzyme concentration is etotmt= m. results of the wild-type estimate appear in the previous section. additional data for the mutant were generated as above based on the c parameter values of mutant/drug treated sample and simulated data sets are used. the estimates for mutant/drug treated sample are shown in table c parameter estimates for acetylornithine aminotransferase from mutant/drug treated sample. results are based on simulated data sets c 1 c 2 c 3 to obtain the prediction of total enzyme change, we take cimt/ciwt. the results are shown in table e tot change prediction based on simulated data sets etotwtetotmt=c1mtc1wt etotwtetotmt=c2mtc2wt etotwtetotmt=c3mtc3wt since any of the ci parameters can be used to estimate the change in etot, the one that gives minimum standard error, c was chosen. this also demonstrates that even though sometimes identifiability issues can occur and some parameters cannot be estimated, our method could still be very useful if one parameter can be estimated accurately. discussion this work is intended to be a first step towards estimating parameters for reactions in large metabolic networks in vivo. in vivo estimation will need to be based on data that have relatively large relative errors in all measured parameters, and will have to deal with a variety of reaction kinetics, including reactions that are reversible and have multiple substrates and/or products. although measurement and estimation of enzyme kinetics has been studied for many decades, there is no single existing estimation method that addresses all of these issues. we have proposed a maximum likelihood approach to estimate kinetic parameters using nonlinear optimization, with estimates on the standard error and bias of the results using the bootstrap. simulations show that invest produces accurate estimates for realistic high-throughput metabolomics data. for example, with 20– samples with coefficients of % in metabolite concentrations and % in reaction rate estimates, estimates have a relative standard error of less than %. collecting data of this quality would be technically difficult, but is within the current state of the art. an advantage of the method is that it estimates each set of reaction parameters independently. if measurements are not available for some metabolites, it can still estimate parameters for those reactions for which the data include all substrates and products. solving the problem of in vivo parameter estimation in its full generality will require meeting a number of additional challenges. some reactions have more complex kinetics than those we consider, especially various kinds of inhibition. when the inhibiting metabolite and mechanism of inhibition are known, the approach described here can probably be generalized to accommodate the inhibition mechanism in our future work. otherwise, a process of model selection may be necessary, where competing models are estimated and the quality of the results compared, with appropriate adjustments for model complexity. in addition, it will be necessary to deal with the kinetics of transport reactions, and to take account of different compartments in the cell. parameter identifiability is a difficult issue in in vivo estimation. we have shown that accurate estimates of all parameters require data that is well-distributed over the kinetics curve, but such data will not often be obtainable for several reasons. experimental data must be obtained by perturbing metabolites and fluxes, for example, by adjusting nutrient media, testing mutants, and targeting reactions with drugs. first, accurate estimation may require non-physiological concentrations of metabolites – estimating c for a reaction that is nearly irreversible being an example. more generally, there is usually inadequate controllability of metabolite concentrations and reaction fluxes to obtain the experimental values needed for accurate estimation, for many reasons including concentrations are toxic or inadequate to sustain life, and rate-limiting reactions that make high fluxes in other reactions impossible to obtain. since we can’t estimate everything accurately, it is important to produce estimates of the standard errors of parameter estimates, so we can tell which ones are meaningful. also, as we note above, if some but not all parameters of a reaction can be estimated accurately, the results still may be useful. for example, it is possible to estimate the total concentration or relative activity of an enzyme in wild-type vs. mutant cells when only one of the kinetic parameters is accurately estimated.  in conclusion, a new method, invest, is developed for estimating reaction kinetic parameters in metabolic networks that addresses the specific challenges of in vivo data. invest uses maximum likelihood estimation based on models where all measurements have potentially relative errors. it can be applied to multiple substrate/product reversible enzymatic reactions with a generalized single displacement mechanism. because it is not always possible to obtain good data covering full range of possible metabolite concentrations and reaction rates, certain parameters may be non-identifiable. invest uses bootstrap to estimate the standard errors of parameter estimations that can tell which estimates are reliable. invest enables the estimation of reaction rates directly from concentration measurements. also, comparing the estimated kinetic parameters of a reaction in the wild type and mutant cells can quantify enzyme abundance or activity change due to genetic variants. the same method can also be used to measure the effect of a drug that targets a particular enzyme. moreover, estimated parameters can be used to estimate maximum reaction rates, which could be used as constraints to improve flux-balance analysis. additional files additional file experimental data support on relative error model. figure s noise errors in high-throughput metabolomic data tend to be relative. the plot shows the empirical standard deviation vs. mean of metabolite concentrations in a publicly available mass spectrometry data set of human urine samples. each sample has technical replicates, which were used to calculate the standard deviation and mean of metabolite concentrations. the data for “peak 105” were chosen because the chromatographic peak appears in all three replicates of the sample and the measurements cover a wide range of concentrations across different samples. low concentrations are omitted because they are highly inaccurate due to background noise. there is a linear relationship between standard deviation and concentration mean, showing that errors are proportional to measured concentration. additional file derivation for steady state rate equation. this file provides a detailed derivation for steady state rate equation of a single reactant and single product reversible metabolic reaction presented in methods section. additional file equilibrium constant k eq. the equilibrium constant k eq is assumed to be a known constant. this file provides the calculation of equilibrium constant k eq based on standard gibbs free energies of formation. additional file standard gibbs free energy of formation metacyc. this file provides standard gibbs free energy of formation taken from metacyc database for metabolites used in the simulation example in results section. availability of data and materials matlab code for reproducing the results of this work is available at stanford digital repository, http://purl.stanford.edu/bg158sn competing interests the authors declare that they have no competing interests. author’s contributions wz and dld defined the problem. wz did the mathematical derivations, implemented the method, and performed the simulations. rk suggested the maximum likelihood approach and worked out key ideas about the approach. dld oversaw the work. the authors collaborated on the writing and all authors have read and approved the final manuscript. 
invest uses the bootstrap method to evaluate the accuracy of its estimates. a new method is needed to estimate kinetic parameters taking into account these factors. a new method, invest, is described for estimating reaction kinetic parameters, which addresses the specific challenges of in vivo data. this data could potentially be used to estimate kinetic parameters for many metabolic reactions. however, high-throughput in vivo measurements have special properties that are not taken into account in existing methods for estimating kinetic parameters, including significant relative errors in measurements of metabolite concentrations and reaction rates, and reactions with multiple substrates and products, which are sometimes reversible. high-throughput assays such as mass spectrometry have opened up the possibility for large-scale in vivo measurements of the metabolome. it has potential applications of quantifying the effects of genetic variants, inference of the target of a mutation or drug treatment and improving flux estimation. simulations show that invest produces accurate estimates for a reversible enzymatic reaction with multiple reactants and products, that estimated parameters can be used to predict the effects of genetic variants, and that invest is more accurate than general least squares and graphic methods on data with relative errors. invest uses maximum likelihood estimation based on a model where all measurements have relative errors. invest addresses several challenges of in vivo data, which are not taken into account by existing methods. relative errorenzymatic reactionparameter estimationmaximum likelihooderror-in-all-measurementsin vivo dataissue-copyright-statement© the author 2015 when data have relative errors, invest produces more accurate and robust estimates.
 high-throughput assays such as mass spectrometry have opened up the possibility for large-scale in vivo measurements of the metabolome. this data could potentially be used to estimate kinetic parameters for many metabolic reactions. however, high-throughput in vivo measurements have special properties that are not taken into account in existing methods for estimating kinetic parameters, including significant relative errors in measurements of metabolite concentrations and reaction rates, and reactions with multiple substrates and products, which are sometimes reversible. a new method is needed to estimate kinetic parameters taking into account these factors.  a new method, invest, is described for estimating reaction kinetic parameters, which addresses the specific challenges of in vivo data. invest uses maximum likelihood estimation based on a model where all measurements have relative errors. simulations show that invest produces accurate estimates for a reversible enzymatic reaction with multiple reactants and products, that estimated parameters can be used to predict the effects of genetic variants, and that invest is more accurate than general least squares and graphic methods on data with relative errors. invest uses the bootstrap method to evaluate the accuracy of its estimates.  invest addresses several challenges of in vivo data, which are not taken into account by existing methods. when data have relative errors, invest produces more accurate and robust estimates. invest also provides useful information about estimation accuracy using bootstrapping. it has potential applications of quantifying the effects of genetic variants, inference of the target of a mutation or drug treatment and improving flux estimation. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords relative errorenzymatic reactionparameter estimationmaximum likelihooderror-in-all-measurementsin vivo dataissue-copyright-statement© the author 2015 