metabolic pathway models consist of enzymatic reactions described by their stoichiometry , the enzymatic rate laws , and their kinetic constants ( such as , for instance , equilibrium constants or catalytic constants ) .
the more we know about these quantities , the more reliably we can simulate the metabolic dynamics .
kinetic laws of individual enzymes have been studied experimentally for about 100 years,(1 ) and metabolic control theory,(2 ) a theoretical apparatus for the analysis of metabolic systems , has been developed since the 1970s .
recently , comprehensive web databases , advances in high - throughput experiments , and inexpensive computing power have led to a new interest in metabolic modeling .
in particular , the numerous large - scale metabolic networks reconstructed from sequenced genomes call for automatic routines that can fill these networks with enzymatic rate laws and turn them into dynamic models .
unfortunately , the enzymatic mechanisms and the rate laws of most enzymes are unknown , and it is laborious to determine them exclusively by enzyme assays .
a pragmatic solution is to substitute missing kinetic laws by standard rate laws , such as mass - action kinetics , generalized mass - action kinetics,(6 ) or linlog kinetics . here , we will use the common modular rate law,(9 ) a generalized version of the reversible michaelismenten rate law , suitable for any reaction stoichiometry and accounting for various types of allosteric regulation .
once a metabolic network and enzymatic rate laws have been chosen , we need numerical values for the kinetic constants . this can be a challenge , especially for large networks .
modelers can find known kinetic constants in published models , in the literature , or in public web resources such as sabio - rk,(10 ) brenda,(11 ) and nist . as pointed out by alberty,(14 ) varying conditions such as ph or salt concentrations
can be taken into account by describing biochemical reactants and reactions in terms of transformed thermodynamic quantities . in the future , automated enzyme assays might provide more kinetic data , but they will still not reach the speed at which metabolic networks are reconstructed from newly sequenced genomes .
available kinetic data may not be suited for a model if they are contradictory or measured under inappropriate conditions ( e.g. , ph values and temperatures ) .
furthermore , data collected from various sources are very unlikely to represent a thermodynamically consistent set .
since incompleteness of the kinetic constants remains a major obstacle , methods for guessing unknown kinetic constants or adjusting the known values will become increasingly important . here , we present parameter balancing , an approach to infer complete and consistent sets of model parameters from incomplete , inconsistent kinetic data .
this is only possible due to mutual dependencies between the kinetic constants and other model parameters , arising from their definitions or from thermodynamic laws ( wegscheider conditions(15 ) and haldane relationships ) . in a simple approach ,
incomplete kinetic data sets could be complemented by inserting all available values into the model and adding other quantities that can be directly computed from them . however , this might leave parameters undetermined and would not eliminate inconsistencies between the original data values . as a better strategy , we determine parameter sets that are consistent by construction and resemble the original data as closely as possible . since these values may not be uniquely determined , we have to restrict them to plausible ranges and quantify the remaining uncertainties by error bars .
we have previously suggested(16 ) implementing this parameter balancing approach in a bayesian statistical framework(17 ) and based on standard rate laws.(9 ) the critical step is to identify a set of independent model quantities that can be freely selected during parameter fitting , sampling , or optimization and from which all remaining quantities can be easily computed . after establishing this dependence scheme for a wide range of kinetic constants and metabolic quantities , we obtain a relatively general and simple data integration method that can guarantee consistent parameter sets . here
, we use more general formulas and present an interactive online workflow for models in sbml format ( systems biology markup language,(18)www.sbml.org ) , then illustrate its use with an example case .
parameter balancing exploits the fact that kinetic model parameters often share mutual dependencies , either by their definition or because of thermodynamic constraints.(16 ) for kinetic models with modular rate laws(9 ) and many other rate laws , the model parameters can be split into two subsets : a set of independent basis quantities that can be chosen arbitrarily and a set of derived quantities that can be computed from linear combinations of the basis quantities . as a simple example ,
let us consider the concentration and the amount of a substance within a cellular compartment .
the amount a can be computed from the concentration c and the compartment volume v by the formula a = cv . if we choose ( e.g. , guess , sample , fit , or optimize ) all three model parameters independently , this dependence will be violated .
of course , this problem is easy to avoid : we just need to treat a as a derived quantity to be computed from the basis quantities c and v. on a logarithmic scale , we obtain the simple linear dependence scheme ln a = ln c + ln v. in parameter balancing , we do exactly the same thing , but treat all quantities of a ( possibly large ) kinetic model at the same time .
in addition , we consider not only dependencies arising from quantity definitions , as in this simple example , but also more involved dependencies arising from the laws of thermodynamics .
the resulting dependence scheme consists of many linear equations , emerging from a thermodynamic and kinetic analysis of the rate laws . for practical calculations ,
these equations are represented by a large , sparse matrix to be constructed from the metabolic network .
initially , the basis quantities are estimated from data by a linear regression model that follows directly from the dependence scheme . following this ,
the estimation is based on bayesian statistics , combining the data with prior distributions for all model parameters , which can be selected to incorporate general knowledge about plausible values .
accordingly , parameter balancing yields not only point estimates but also posterior distributions , from which mean values , variances , correlations , and even random samples of all relevant quantities can be obtained .
on the basis of our practical experience , we have extended the original parameter balancing approach in three major ways : ( 1 ) our model quantities comprise not only kinetic constants , but possibly also metabolite and enzyme concentrations for one or more metabolic states . from their values , we can derive other state - dependent quantities ; in particular , the chemical potentials and reaction affinities .
the reaction affinities describe how strongly chemical reactions deviate from their equilibrium and predetermine the reaction directions .
( 2 ) in biochemical reactions , individual protonation states of a molecule (  chemical species  ) are usually subsumed in a single  reactant  concentration .
as pointed out by alberty , these reactants should be described by transformed thermodynamic quantities , which effectively depend on the ph value and on salt concentrations . in our approach , quantities such as equilibrium constants , gibbs free energies , and chemical potentials are given as transformed values , depending both on temperature and on ph . if experimental values stem from measurements under incompatible conditions , the discrepancies can be corrected during parameter balancing .
( 3 ) although prior distributions can help to define plausible ranges for the basis quantities , they are not applicable to the derived quantities . as a consequence ,
whenever few data on these quantities are available , their large uncertainties lead to unreasonable balancing results .
we address this problem by augmenting the experimental data set with fictitious pseudo values , playing a role similar to the prior distributions .
they allow the modeler to control the variance of the independent values and , thus , the reliability of the whole estimate . a detailed description of the original method and all new features is given in the supporting information .
for practical use , we have implemented an interactive workflow for sbml models , allowing the user to balance the parameters and to replace or complete kinetic rate laws in the model .
three kinds of information are needed as an input : the network structure , which is obtained from the sbml file , the mathematical rate laws , which are chosen from the list of modular rate laws,(9 ) and a table with collected kinetic constants and other relevant data ( for an example , see figure 1 and table 1 ) . in the model , enzymes and allosteric activators and inhibitors
the quantity types are defined as in the systems biology ontology,(21 ) and the names of the reactions and species refer to ids in the sbml model .
all this helps the process to run in an almost fully automated manner . for further details ,
the model structure is defined by the sum formula f6p + atp  fbp + adp with the molecular species f6p ( fructose 6-phosphate ) , fbp ( fructose 1,6-bisphosphate ) , atp ( adenosine triphosphate ) , and adp ( adenosine diphosphate ) . in our parameter balancing workflow ,
the stoichiometry is provided as an sbml file , and kinetic constants ( shown in flags ) are read from a separate data file ( see table 1 ) .
the constants shown suffice to define a common modular rate law,(9 ) which we use as a standard rate law . at the end of the workflow , kinetic rate laws with a set of balanced parameters are inserted into the sbml model . as a result of balancing , most parameters are represented by normal distributions of their logarithmic values . when converting these values back to nonlogarithmic scale , we obtain log - normal distributions and need to carefully distinguish between their arithmetic and geometric mean values . to insert the most probable parameter set into a model , we choose the median values of the nonlogarithmic parameters , which are identical to the geometric mean .
in addition , we can sample additional logarithmic parameter sets from the normal distribution , convert them back by applying the exponential function , and insert them into the model . in the supporting information
, we discuss a number of possible extensions , such as handling of identical species in different compartments , electrochemical potentials , cell compartments with different ph values , forward and reverse reaction rates , use of correlated priors , use of equilibrium constants as basis quantities , and prescribed reaction directions .
the web site also contains the documentation and a number of example models , including the phosphofructokinase reaction discussed below . at the beginning of the workflow
, the user uploads an sbml model,(18 ) defining the network structure ( see figure 1 ) as well as a formatted data table ( see table 1 ) listing the known kinetic constants and other quantities relevant for the model .
after uploading both files , the data can be filtered for a specific source organism , and missing standard errors are completed by default values .
then , a table of relevant model quantities is produced , with blank rows where data are missing and averaged values where more than one data point was available .
the standard chemical potentials or equilibrium constants measured under different temperature or ph are not averaged , but kept as separate values . as a second source of information , we use prior distributions , describing our general expectation about the quantity types .
such priors can be derived from the typical ranges of kinetic constants found in the literature . on the basis of the prior distributions and pseudo values chosen by the user , a set of model parameters
, the user can also choose between several simpler methods for data completion : ( i ) completing all missing quantities by the prior means and widths ( for missing basis quantities ) or by pseudo values and their standard errors ( for missing derived quantities ) , leading to complete , but inconsistent parameter sets ; ( ii ) completing all missing basis quantities by prior values and computing all derived quantities by the dependence scheme ; ( iii ) completing all missing derived quantities by pseudo values , followed by balancing without pseudo values ; or ( iv ) balancing without pseudo values .
afterward , the user can choose a rate law from the list of modular rate laws(9 ) that will be inserted into the sbml file .
the parameter values represent either median values or random samples from the posterior distribution . finally , the balanced quantities and the completed model are exported , respectively , as a data table and as a fully parametrized sbml file .
we have tested parameter balancing with medium - scale models of central metabolism : yeast glycolysis models by teusink et al.(22 ) and hynne et al.(23 ) and a model of metabolism in pancreatic beta cells.(24 ) the original models , the collected data , and the balancing results for all models can be found at www.semanticsbml.org in the online documentation . for simplicity , we consider here a small model of the phosphofructokinase reaction , a key step in upper glycolysis : 
  the enzyme phosphofructokinase ( pfk ) , which transfers a phosphate group from atp ( adenosine triphosphate ) to fructose 6-phosphate ( see figure 1 ) , has been studied extensively , but the databases brenda,(11 ) nist,(13 ) and sabio - rk(10 ) do not contain a complete kinetic data set for a reversible rate law .
therefore , we pretend that the kinetic law of pfk is unknown , apply parameter balancing to the pooled data from several organisms , and compare the results to parameters from published kinetic models for the baker s yeast saccharomyces cerevisiae . our sbml model for the pfk reaction , including miriam - compliant annotations , was automatically constructed from the kegg(27 ) reaction identifier r04779 with the tool semanticsbml(28 ) ( accessible at www.semanticsbml.org ) .
then , we collected data from the databases sabio - rk,(10 ) brenda,(11 ) nist,(13 ) and yeastgfp(29 ) and from publications by nissler et al.,(30 ) albe et al.,(31 ) and alberty.(20 ) the preprocessed data are shown in table 1 . for the concentrations of fructose 6-phosphate and fructose 1,6-bisphosphate , we inserted pseudo values ( arithmetic mean values arising from geometric mean values 0.5 and 1 mm , respectively , with a broad standard deviation of 0.5 for the decadic logarithms ) .
some of the values were obtained by averaging over several data values measured in different organisms .
the set of balanced parameters , obtained with the default workflow settings , is shown in table 2 . for abbreviations , see figure 1 .
data taken from alberty,(20 ) yeastgfp,(29 ) nissler et al.,(30 ) nist,(13 ) brenda(11 ) . a comparison with values from the literature and existing models
the maximal velocity of the phosphofructokinase , mainly determined by a broad prior on catalytic constants and by an experimentally derived count number of the enzyme molecules , remains within a sensible range of the values from the literature ( 0.006 mm / s ) .
we find the equilibrium constant to be much lower ( 0.072 ) than the one estimated by teusink et al.(22)(80 ) , which is clearly due to the small input value ( 0.08 ) from nissler et al.(30 ) finally , the balanced inhibition constant for atp ( 0.376 mm ) lies within a sensible range of the one found in sabio - rk(10 ) ( 1.09 mm ) , but is nonetheless a little lower , since the input data value ( 0.396 mm ) is lower , as well .
the equilibrium constants play a central role in parameter balancing and their numerical values primarily arise from measured values and from known standard chemical potentials . to incorporate their dependence on ph values , we reran the balancing with input data containing equilibrium constants measured at different temperatures and ph values ( see table 5 ) .
instead of just averaging over them ( as for other duplicate values ) , parameter balancing can adjust the values to a certain target temperature and ph value chosen by the user .
the aim is , of course , to describe in vivo conditions considered in the model . when choosing a target temperature of 300 k and a target ph value of 7 , we obtained a balanced equilibrium constant k = 0.19397 . since the input measurement conditions ( average ph , 7.7 ; average temperature , 303.82 k ) differ from the desired conditions ( ph , 7 ; temperature , 300 k ) , the input value for the equilibrium constant ( k = 0.02923 ) changes significantly after parameter balancing .
when we choose target conditions closer to the data measurement conditions ( ph , 7.7 ; temperature , 304 k ) , the balanced value k = 0.19393 moves slightly closer to the data value . finally , we tested what would happen without any direct information on equilibrium constants .
when we remove all data and pseudo values for the equilibrium constant , we obtain abnormally high values for the majority of the derived parameters .
even though equilibrium constants can in principle be obtained from standard chemical potentials , this result indicates that the remaining uncertainty may be extremely large , and pseudo values are an efficient way to limit the ranges of balanced values .
parameter balancing can also be applied to models of bigger size . due to the large amounts of data that are produced , the result tables are not included in this publication . instead , the detailed kinetic data for the models of teusink et al.(22 ) ( 17 reactions , computing time 0.29 s ) , hynne et al.(23 ) ( 24 reactions , computing time 0.50 s ) , and jiang et al.(24 ) ( 45 reactions , computing time 2.32 s ) can be found , downloaded , and used on our web site www.semanticsbml.org .
nevertheless , the rising resource demands of bigger models can be a limiting factor . in the supporting information
the example of the phosphofructokinase reaction shows that parameter balancing can produce consistent estimates on kinetic constants in plausible ranges and close to available input data . in the spirit of bayesian statistics
, we do not estimate the unknowns by averaging over the data , but by fitting the data with a statistical model ; in our case , produced from the dependence scheme .
one advantage is that such a model can handle not only uncertainties of individual quantities , but also correlated uncertainties arising from parameter dependencies .
of course , this approach strongly depends on the collected input data and on the chosen prior distributions . in situations when input data are missing
as a side effect , balancing will shift all the values  even if data are available for a certain quantity  toward the center of the prior distribution and toward the pseudo values , as can be seen in the example shown in tables 3 and 4 .
if this effect is unwanted , there are two possibilities to avoid it : on one hand , experimental values can be fixed by assigning small standard errors to them ; on the other hand , our workflow also provides the possibility to replace unknown values by prior or pseudo values without further balancing .
the median posterior values obtained from balancing can be used as point estimates , but we can also sample parameter sets from the posterior distribution . studying models with sampled parameters can be an emergency solution if few data are available , but it is also a convenient way to explore the dynamics in metabolic networks for a wide range of kinetic parameters and to discern the influences of network structure and kinetics on the dynamic behavior .
parameter balancing is a general approach that could be applied to any physical model as long as all parameters fit into a linear dependence scheme . for kinetic models
, we could establish linear dependencies for most relevant quantities , considering some of them on logarithmic scale .
the choice of model parameters to be balanced is not fixed , but can depend on the specific situation ( for a variety of possible variants , see the supporting information ) . in the scheme presented here ,
the equilibrium constants are derived from standard chemical potentials , which allows us to incorporate data or predictions of gibbs free energies of formation ; a general alternative , with even fewer parameters , would be to express all equilibrium constants by a set of independent equilibrium constants ( see the supporting information ) . however , since model identifiability is guaranteed by the usage of prior distributions , keeping the number of parameters small is usually not crucial .
some of the constants ( for instance , the inhibition constants ) are independent of all other parameters and could therefore be balanced separately , but as long as the parameter set is not too large , it turned out to be practical to account for all kinetic constants and metabolic data in one large dependence matrix .
the reaction rates would be the most important target but do not fit into the scheme . in the future , their signs and the individual forward and
once a subset of the flux directions have been predefined , available data on metabolite levels can directly contribute to the balancing of kinetic constants and thereby improve the estimation results .
arguably , the most critical point in our approach is the use of heterogeneous data from different sources .
ideally , all kinetic constants used in a model should stem from measurements under the same , standardized , nearly in vivo conditions .
since such data are rarely available , modelers often need to utilize heterogeneous data collected from literature and databases , acknowledging that this may cause various problems . first , kinetic constants measured in vitro and in vivo differ due to different ph values , temperatures , salt concentrations , or other factors .
although we attempt to take these conditions into account , it may be difficult to apply corrections , since the exact conditions in living cells are not known .
second , most kinetic models neglect the complexity of the cell ( e.g. , molecular crowding , channeling , inhomogeneous localization of enzymes ) .
since the kinetic parameters in such models describe an effective behavior ( e.g. , averaged over different cell compartments ) , they will differ from the values measured in vitro .
third , there are different conventions for reaction formulas ( e.g. , h ions and h2o molecules may or may not be listed ) and about the standard concentration c ( in our case , 1 mm ) . since the definition of transformed equilibrium constants crucially depends on the reaction formulas , it is important that model and data sources use the same , appropriate conventions . finally ,
if kinetic constants are taken from existing models , their values may become invalid out of this specific context .
thus , especially for poorly identifiable parameters , it is important to consider the uncertainties in the original estimations .
facing these difficulties , parameter balancing follows a pragmatic approach : even if data have a low quality , we may still use them as clues about unknown parameters .
if a quantity has been measured under the wrong conditions , we may account for this by increasing its standard error .
on one hand , this will decrease the weight of this data point in the balancing process ; on the other hand , the uncertainty of all related balanced parameters will increase , reflecting our precautionary approach . in some cases , we may need to assign very small uncertainties to some of the input data .
for instance , if an sbml model contains kinetic laws and we intend to replace only some of them , we have to make sure that the new rate laws are compatible with those pre - existing rate laws that will remain in the model . as a key precondition , all equilibrium constants in the resulting model need to satisfy the wegscheider conditions.(15 ) to ensure this , one can determine the equilibrium constants of the existing rate laws and use them as input data ( with zero standard error ) in parameter balancing .
due to the large uncertainties in standard chemical potentials and metabolite concentrations , it is unlikely that models obtained from parameter balancing will directly show realistic stationary flux distributions . in the future , this could be enforced by a subsequent fit to  omics 
data(16 ) or by predefining the signs of reaction affinities . these signs will induce dependencies between kinetic and metabolic quantities ; for instance , between the standard chemical potentials and the metabolite concentrations . if several metabolic states are treated within the same dependence scheme , the resulting method would resemble network - embedded thermodynamic analysis(32 ) and allow to use the results of previous fluxome and metabolome analysis ( by flux - balance methods including thermodynamic constraints ) as input data for parameter balancing .
parameter balancing yields consistent and complete parameter sets for kinetic models , as a potential starting point for further modeling .
it can integrate incomplete and contradictory input data and respects constraints implied by common modeling assumptions ( standard rate laws ; reactants in ideal , dilute solution ) .
if few input data are available , parameter balancing can help to quantify the remaining uncertainties , whereas with more and higher - quality data , the predictions become more trustworthy and precise .
the resulting posterior distribution can be used to define parameter ranges , to sample possible parameter sets , or to be reused as a prior for following rounds of parameter balancing .
 kinetic modeling of metabolic pathways has become a major field of systems biology . 
 it combines structural information about metabolic pathways with quantitative enzymatic rate laws . 
 some of the kinetic constants needed for a model could be collected from ever - growing literature and public web resources , but they are often incomplete , incompatible , or simply not available . 
 we address this lack of information by parameter balancing , a method to complete given sets of kinetic constants . 
 based on bayesian parameter estimation , it exploits the thermodynamic dependencies among different biochemical quantities to guess realistic model parameters from available kinetic data . 
 our algorithm accounts for varying measurement conditions in the input data ( ph value and temperature ) . 
 it can process kinetic constants and state - dependent quantities such as metabolite concentrations or chemical potentials , and uses prior distributions and data augmentation to keep the estimated quantities within plausible ranges 
 . an online service and free software for parameter balancing with models provided in sbml format ( systems biology markup language ) is accessible at www.semanticsbml.org . 
 we demonstrate its practical use with a small model of the phosphofructokinase reaction and discuss its possible applications and limitations . in the future 
 , parameter balancing could become an important routine step in the kinetic modeling of large metabolic networks . 
Introduction
Methods
Results
Discussion
Conclusions