BACKGROUND
metabolic activities of living cells are accomplished by a well regulated, highly coupled network of numerous enzyme catalyzed reactions and selective membrane transport systems.
to engineer such systems, enzymatic, transport and regulatory functions of the cells are manipulated via the use of recombinant dna technology  <cit> .
within the purpose of metabolic engineering, i.e.
rational redesign of the metabolic systems, the highly relevant question of which  perturbation should be applied in order to increase the productivity of the microorganism is addressed.
the answer to this question requires information on both the regulatory level and the metabolic level.
ter kuile and westerhoff showed by their hierarchical control analysis that the pathway flux is rarely controlled solely by gene expression, but that metabolite levels are also relevant  <cit> .
in this context, we focus on the metabolic reaction network level.
modeling metabolic reaction systems is usually based on stoichiometric sometimes followed by kinetic modeling.
one of the initial steps in the modeling of metabolic reaction networks is to determine the structure and steady state characteristics of a given network using stoichiometric information alone.
steady state models describe time invariant fluxes, gathered from steady state experiments; hence they reflect the structural characteristics of the system.
metabolic flux analysis  and metabolic network analysis  were developed as powerful tools to analyze such flux data.
at steady state, the mass balances over the metabolites in the metabolic network yield a set of linear relations between the metabolic fluxes which can be expressed as:
s·v =  <dig>
where s is the  stoichiometric matrix and v is the  vector of metabolic fluxes, where m is the number of balanced metabolites and r is the number of fluxes.
here, the system is in most of the cases underdetermined so that there are an infinite number of possible solutions.
the realized solution depends on the kinetic properties of the involved reactions; this information is seldom known.
to bypass the need of information on kinetics of individual reactions, alternative mathematical approaches have been proposed in the past to obtain a unique solution.
an example is the constraint based optimization approach which is based on assumed optimality criteria, e.g.
maximum growth, given biochemical, thermodynamic and irreversibility constraints and maximal reaction rates  <cit> .
later, segre et al.
proposed the optimality constraint that requires maximization of biomass formation while minimization of metabolic adjustment  in order to obtain a unique flux distribution of a mutant strain.
in their approach, they defend their optimality criterion that a knock-out mutant strain would optimize its biomass production rate by changing minimally its metabolic fluxes from the wild type strain  <cit> .
an alternative metabolic modeling framework which uses a fitness function is the cybernetic approach.
this approach assumes that an organism is an optimal strategist in utilizing all available sources with maximum efficiency.
the expression and activity of the enzymes that catalyze network functionality are regulated by cybernetic control variables obtained from the solution of a constrained optimization problem  <cit> .
despite a number of successful applications especially in mixed substrate and prediction of knock-out lethality, all stoichiometric modeling approaches have their limitations, e.g.
they can not predict time courses of the cellular processes, are based on "assumed goals" of the cell, and do not give insight in molecular events occurring in the cells, since information about the kinetic properties of the individual enzymes are not required.
in the light of the above arguments, it is apparent that, to advance our understanding of the functioning of living cells a systems biology approach is needed, whereby the use of dynamic mathematical models of metabolic reaction networks to describe the complex kinetic behavior and interactions  is becoming increasingly relevant.
the kinetic behavior of many important enzymes occurring in metabolic networks have been studied extensively, however these studies have generally been performed under non-physiological conditions in test tubes , and therefore the applicability of these results to the in vivo metabolism is doubtful  <cit> .
teusink et al.
showed that discrepancies exist between the in vivo measured changes of the concentrations of the glycolytic metabolites and their estimates using models based on mechanistic rate equations and in vitro determined parameters  <cit> .
this basic problem invalidates detailed models of metabolism containing kinetic parameters which have been determined in vitro.
therefore, it is preferred to base the kinetic analysis of metabolic networks on in vivo studies of intact cellular networks.
these in vivo studies are based on steady state and/or dynamic perturbations of a metabolic network starting from a reference steady state that is defined by its fluxes, enzyme activities, metabolite levels, extracellular concentrations, and kinetic parameters.
it is also important to notice that the number of parameters that is typically involved in the traditional mechanistic equations is very large which causes identifiability problems due to parameter insensitivity.
despite the information richness of data obtained from dynamic perturbations in in vivo experiments, there is a limit for identification of the parameters.
for these reasons, it seems justified not to limit ourselves to the available complex mechanistic enzyme kinetics which must also be considered as approximations of the true in vivo behavior  <cit> .
approximative enzyme kinetic formats, which contain much fewer kinetic parameters are therefore of general interest.
an overview of different approximative kinetic formats  used in metabolic network modeling is given by heijnen  <cit> .
one of the proposed formats is, linlog kinetics, which has been introduced for modeling of in vivo kinetics and for metabolic redesign, and shown to have a good approximation quality, standardized format and relatively few parameters  <cit> .
in linlog kinetics, all the rate equations have the same mathematical structure in which the relation between rates and enzyme levels is proportional, while for metabolite levels, a linear sum of logarithmic concentration terms is proposed.
all variables are considered relative to a reference steady state ).
the linlog approximation is valid in the neighborhood of the reference state ), but quite large changes of metabolite concentrations, enzyme levels and fluxes are allowed  <cit> .
the parameters  in the kinetic equations are the same scaled elasticities
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaqadaqaaggaciab=v7alnaadaaaleaacqwg4baeaeaacqwg2bgdaagccqgh9aqpdawcaaqaaiabdiha4naacaaaleqabagaegimaadaaagcbagaemosao0aawbaasqabeaacqaiwaamaaaaaowaaqgaaeaadawcaaqaaiabgkgi2kabdaha2bqaaiabgkgi2kabdiha4baaaiaawia7amaabaaaleaacqwg4baecqgh9aqpcqwg4baedaahaaadbeqaaiabicdawaaaasqabaaakiaawicacaglpaaaaaa@45b3@ that are used in metabolic control analysis .
it is important to note that the elasticity parameters appear in the model in a linear fashion.
vj0=ee0+εcvln⁡)
mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaiabdaha2bqaaiabdqeaknaacaaaleqabagaegimaadaaaaakiabg2da9maalaaabagaemyzaugabagaemyzau2aawbaasqabeaacqaiwaamaaaaaowaaewaaeaacqaixaqmcqghrawkiigacqwf1oqzdaqhaawcbagaemieaghabagaemodayhaaogagiibawmaeioba42aaewaaeaadawcaaqaaiabdiha4bqaaiabdiha4naacaaaleqabagaegimaadaaaaaaogaayjkaiaawmcaaiabgucariab=v7alnaadaaaleaacqwgjbwyaeaacqwg2bgdaagccyggsbabcqggubgbdaqadaqaamaalaaabagaem4yamgabagaem4yam2aawbaasqabeaacqaiwaamaaaaaagccagloagaayzkaaaacagloagaayzkaagaaczcaiaaxmaadaqadaqaaiabikdayagaayjkaiaawmcaaaaa@5758@
when the elasticities are known, a full dynamic model of the whole metabolic network can be set up using linlog kinetics.
such a model allows in principle the calculation of control coefficients also under dynamic conditions.
in linlog kinetics, the elasticities are the kinetic parameters represented in the elasticity matrix.
from these, and a given network structure specified in the form of a stoichiometric matrix, the control coefficients for a reference condition can be calculated from the summation and connectivity relations developed in the framework of mca.
also the change in control coefficients upon large changes in enzyme levels can be calculated  <cit> .
moreover, the linlog formulation enables the analytical solution of the mass balances for steady state metabolite and flux levels in the metabolic network, providing the solution of the metabolic redesign problem, i.e.
determination of the optimal enzyme levels that maximize a certain flux while the total amount of enzyme and the metabolite levels are constrained.
visser et al., reported a successful application of linlog kinetics in an in silico study that aims to solve this metabolic redesign problem  <cit> .
in order to determine the kinetic parameters of a model of a given in vivo metabolic system, in vivo perturbations of the complete metabolic network have to be performed.
there are two main types of perturbations that can be imposed on the system: steady state and dynamic perturbations.
in steady state perturbations, usually the enzyme activity of one  of the reactions is changed by adding specific inhibitors or activators or by genetically changing the enzyme activity, resulting in a new steady state.
in steady state perturbations an important problem has been addressed by kacser and burns, which is the determination of the set of reactions that has to be perturbed in steady state fashion, in order to be able to determine all elasticities for a given metabolic network which resulted in their "double modulation"  <cit> .
they showed that for a simple linear chain of reaction, perturbation of the activities of the first and last enzymes allows determination of the elasticities of all enzymes of the chain under the condition that each enzyme is only responsive to its substrate and product, which rules out the fact that feedback loops are present.
the theoretical basis is presented in later studies which showed that determination of the elasticities for any enzyme in such a simple chain requires two perturbations, one upstream and one downstream of the enzyme concerned  <cit> .
giersch and cornish-bowden extended the double perturbation approach to more complex pathways containing branch points, regulatory loops, and conserved moieties and they provided guidelines to list the possible reactions to be modulated in order to determine the elasticities for arbitrary metabolic networks  <cit> .
from the obtained list a minimum set of steady state perturbations, to be imposed on a specific network, can be chosen.
as an alternative to the analysis of giersch and cornish-bowden, hofmeyr and cornish-bowden offered co-response analysis, to identify the mono-functional units that respond together to any perturbation applied  <cit> .
these mono-functional units have to be dissected in order to determine the elasticities belonging to these groups.
linlog kinetics has been successfully applied to estimate the elasticity parameters for a linear pathway from sets of steady state metabolite concentrations and enzyme activities  <cit> .
additionally, heijnen et al.
proposed a method to obtain flux control coefficients around a branch point, from large enzyme perturbation leading to large steady state flux perturbations, using the linlog kinetic format  <cit> .
this approach allows obtaining explicit solutions for steady state flux and metabolite levels as a function of large changes in enzyme activities.
an alternative to steady state perturbations are dynamic perturbations, in which the system, being initially at the reference state, is disturbed to create time dependent data of transient metabolite levels.
the dynamic metabolite profiles are typically obtained as a series of snapshots in time.
a recent example of these are in vivo measurements of a number of metabolites in rather dense time sequences of a few seconds or minutes, using 'rapid sampling' methods with subsequent high precision metabolite measurement techniques  <cit> .
transient data are rich in information and allow determination of the time hierarchy of the different elements of the metabolic network and the causal relationships between the network elements.
it is possible to exploit these data to estimate the parameters of a traditional full kinetic model  <cit> .
such a model can subsequently be used to calculate the values of the elasticities at a given reference steady state.
alternatively, when the linlog kinetic format is used the elasticities can be directly estimated from the transient data.
this was demonstrated by kresnowati et al.
who estimated the elasticities of a small example network from transient concentration data, assuming that the dynamic fluxes are unknown  <cit> .
an important difference between steady state and dynamic perturbations is that in the former both the steady state enzyme levels  and the metabolite concentrations  have to be measured, whereas in rapid dynamic perturbation experiments only xi is required because the enzyme activities can be considered constant within a sufficiently short time window following the perturbation.
note that flux data are required in both methods.
however, fluxes are not independent variables as they follow from the measured intra- and extracellular metabolite concentrations and the proper mass balances.
in this work, we present a method to estimate kinetic parameters , using linlog kinetics, using metabolome data obtained from steady state and dynamic perturbation experiments.
as such, this can be considered as an extension of the work by visser and heijnen  <cit> , wu et al.,  <cit>  and kresnowati et al.,  <cit>  who have provided the theoretical framework and a small practical application of the linlog kinetic format.
we first generalize the notation to make it applicable to networks of arbitrary size and complexity.
additionally, we address further the question of parameter identification from experimental data.
we monitor the propagation of error throughout the proposed parameter estimation procedure and we determine the subset of identifiable parameters.
to illustrate the proof of principle, we apply the presented theory on in silico generated data, with realistic experimental settings, to be able to compare the obtained results with the "known truth".
this data is generated using a previously published glycolytic pathway model  <cit>  of saccharomyces cerevisiae.
RESULTS
in silico pulse experiments
in these experiments, the steady state was perturbed by increasing the extracellular glucose concentration.
when the external glucose concentration is increased, the intracellular glucose concentration also increases as the sugar is transported through facilitated diffusion.
from figure 1a it can be seen that the g6p concentration increases around 6% in the first five seconds and then starts to decrease, as the atp concentration drops due to the phosphorylation of the glucose in the hexokinase reaction.
atp starts to recover after about  <dig> seconds due to increased increase atp production downstream of glycolysis.
this is in agreement with the literature on the glucose uptake dynamics following a sudden increase in the external sugar concentration  <cit> .
following this atp fluctuation, the pep concentration also decreases and then recovers.
in the first five minutes, the concentrations of the external metabolites increase because of their increased production.
after this initial period of  <dig> min, in which the external glucose concentration is nearly constant, the concentration of external glucose slowly decreases due to continuous consumption and wash-out from the reactor.
the cells follow this glucose drop by dropping the levels of intracellular metabolites.
for the extracellular products  the washout from the reactor is larger than their production and hence their concentration also drops.
after about  <dig> minutes, the original steady state is restored .
such a pulse experiment therefore always delivers a highly dynamic dataset, followed by a pseudo steady state dataset.
parameter estimation
the dynamic data obtained during the first five minutes of the glucose pulse experiment  were used to estimate the elasticities via the linlog parameter estimation procedure outlined in the methods section .
in figure 1a, the simulation of the same dynamic perturbation using the estimated elasticities is also given.
in figure  <dig>  the estimated elasticities are compared with the theoretically calculated elasticities, derived from mechanistic rate equations at the reference state.
although the linlog simulation results do not differ much from the noise-free experimental data generated with the mechanistic model, some estimated elasticities are far from the expected theoretical values.
the difference is in some cases in magnitude, for other cases there is a sign contradiction.
specifically,  <dig> out of  <dig> elasticities  are predicted in very good agreement with the theoretical ones ,  <dig> estimated elasticities  agree in sign, but the magnitudes differ by more than 50% , and  <dig> elasticities  are estimated with a sign contradiction.
notice that in the case of e <dig>  the 200% deviation is already expected.
this elasticity explains the effect of g6p on the polysaccharide formation rate and in the mechanistic model this effect is explained with a hill type kinetic equation with a hill coefficient of  <dig> , hence this rate is very sensitive to the changes in g6p.
to mimic this behavior, e <dig> is estimated to be much higher than its theoretical value.
the poor identification of the elasticities has two possible reasons: either the experimental design is poor or there is a structural problem resulting that some interactions in the network can not be resolved from the available information.
from an experimentalist point of view, it can be seen from figure 1a that the changes in the intracellular metabolites during the first  <dig> seconds are between 25–30% which can be easily detected with the current measurement techniques  <cit> .
the extracellular glucose and polysaccharides also change to a detectable extent, but the changes in the ethanol  are only 2–3% which is hard to detect.
to check if the poor experimental design causes problems in the parameter identification, we have altered the experimental design to create larger changes in these concentrations.
first the biomass concentration was increased to  <dig> gdw l- <dig>  also at the same moment that the fermentor was pulsed with the glucose, the inflow and outflow of the fermentor was stopped so the operation was effectively switched from though-flow mode to batch mode, and the extracellular metabolites were not washed out anymore.
this allowed more rapid accumulation of the secreted products resulting in much larger changes.
also, the change in the extracellular glucose concentration was more pronounced because there was no further addition of glucose after the pulse.
these changes did not result in considerable changes in the intracellular metabolite profiles; they follow the slight change in the external glucose profile.
the results of the new experimental design are depicted in figure  <dig> which represents the new experimental data where only the significantly changed external metabolite profiles are depicted, together with the linlog simulation using the estimated elasticities.
however, no significant improvement was achieved in the identification of the elasticities compared to the initial experimental design .
hence it can be concluded that the problem is not due to poor experimental design.
a closer look at the estimated elasticities reveals that the elasticities belonging to the part of the metabolic network, between vhk and vgapd are correctly estimated, whereas the estimation of the elasticities belonging to the lower part of glycolysis and the uptake reaction is poor.
this is mainly due to the fact that the information content of the metabolite data is insufficient to resolve the complex interactions of metabolites and enzymes.
specifically, in the vpk reaction rate, the feed forward activation of fdp and the mass action effect of pep are assumed.
the perturbation of the system via an increase in the external glucose concentration results in an increase in the vpk flux.
however, this does not allow a separate determination as to what extent the change in vpk is due to mass action effect of pep or due to feed forward effect of fdp.
this results in identifiability problems for the elasticity parameters that describe the ethanol production .
the same argument is valid for the elasticities belonging to the glycerol producing reaction, since in the in silico network the glycerol production is assumed to be proportional with the ethanol production; hence the same combined effect is seen in the glycerol production, resulting in an identifiability problem for the latter triplet of elasticities
in order to analyze the information content of the data, the fisher information matrix  is calculated as described in  <cit> .
fim is an indicator of the information content with respect to the parameters and due to the parameter linearity in linlog kinetics fim is independent of the values of the kinetic parameters.
it can be calculated as yty, y being the design matrix appearing in equation .
the singular values of the fim hold information on the number of linear dependencies between the columns of the data matrix.
in our system, we checked the singular values of fim and concluded that not all the singular values are equally significant, which shows that there are some colinearities within the data matrix.
this is also evident from the condition number, which equals the ratio of the highest to the lowest eigenvalues of the matrix and is ideally  <dig> in a completely uncorrelated case.
in this case the condition number was calculated to be  <dig>  ×  <dig>  this fact supports the previous argument that using only one single set of dynamic data is not sufficient to resolve all the elasticity values.
in order to obtain better estimates of the elasticities, steady state perturbation data are introduced.
the steady state enzyme perturbation data are generated as described in the methods section : in order to resolve the feed forward effect of fdp on the vpk reaction, the enzyme amount of the vgapd reaction is perturbed in addition to the vin and vatpase.
using the steady state perturbation data in table  <dig>  and the calculation procedure described in the methods section , we have estimated the elasticities for the in silico network under study.
after the estimation of the elasticities using steady state data only, we have also simulated the dynamic perturbation of figure 1a using linlog kinetics, using the estimated elasticities.
the results are presented in figure  <dig>  figure 4a shows the experimental data and the linlog simulation of the same pulse, and figure 4b compares the elasticities estimated using the steady state data with the theoretical elasticities.
as can be seen from figure 4b, using the steady state data, the elasticities are estimated in good agreement with the theoretical elasticities.
however, as can be observed from figure 4a, the linlog model deviates significantly from the experimental data, so the estimated elasticities have to be further refined.
having a good initial estimate for the linlog elasticities from the steady state perturbation data only, we used in the second step both the dynamic and steady state data for the parameter estimation procedure.
the nonlinear regression procedure allows combining both available steady state and dynamic experimental data, so that all parameters can be accurately estimated.
after the non-linear fit to the data, we have obtained the final estimation of the elasticities.
the simulation of the pulse experiment with these final estimates, and the comparison of the elasticities are given in figure 5a and figure 5b respectively, which shows that most of the elasticities are correctly estimated and that the experimental data are well described.
analysis of parameter identifiability and model reduction
in the previous sections, we have seen that we could not determine all the parameters from dynamic perturbation data only and we have introduced the steady state perturbation data.
after determination of all of the elasticities, using both types of perturbation data, we analyzed back the identifiability of each parameter under noise.
in order to carry the identifiability and error propagation analysis, mc simulation was used.
as described in the methods section , 10% relative error was added to the noise free data represented in figure 1a , and the non-linear estimation procedure was implemented using the elasticities estimated from the steady state data, as the initial guesses.
after repeating this scheme  <dig> times, we obtained a distribution for each of the elasticities from which we calculated the relative error for each elasticity as the standard deviation per mean of the corresponding elasticity .
the adjustable parameter λ in equation  determines the relative importance of the different sources of data.
in our case, these are dynamic and steady state data.
notice that during the mc simulations, noise is added only to dynamic data.
for steady state perturbations, the number of data points that can be obtained is theoretically infinite.
this implies that, considering the central limit theorem, the steady state data can, in principle, be considered as noise free.
however, the aim of the mc simulations was to elucidate the potentially unidentifiable elasticities using dynamic perturbation data, and therefore we chose λ to be equal to zero.
the other extreme  would suppress the effect of noise and wouldn't lead to detection of poorly identifiable elasticities.
the values in between, are up to the choice of the modeler, depending on how many data points have been obtained from a dynamic pulse experiment, the standard deviations of the measurements, etc.
we have implemented different values for the value of λ, but the outcomes of the mc simulation, i.e.
which elasticities can hardly be identified, did not change qualitatively.
from the distribution of each elasticity, the p-value, indicating the probability that the actual elasticity is zero and that the estimated value was caused by the random error only, is calculated via t-test and the results are given in table  <dig> .
the elasticities e <dig>  e <dig>  e <dig>  and e <dig> have a relative standard deviation higher than 200%, so these elasticities are not identifiable under the 10% noise present.
therefore we conclude that these elasticities are clearly candidates for model reduction by assigning them to be zero .
a closer look reveals that e <dig> represents the feed-forward effect of fdp on vpk.
e <dig> and e <dig> represent the effect of atp on vhk and vgol respectively and e <dig> is the feedback inhibition of g6p on vin.
notice that the elasticities of the reactions of vin, vpk and vgol were already badly identified, when we estimated the elasticities using dynamic data only .
using the arguments stated above, we set the values of these four  elasticities to zero, implemented the non-linear regression step, with the new sparser elasticity matrix, and estimated the elasticities using steady state and dynamic data.
the results are presented in figure 6a, where the experimental data and the linlog simulation using the new set of elasticities are given and in figure 6b where the comparison of the estimated elasticities with the theoretical ones is presented.
it is clear that with the new set of  <dig> elasticities  the dynamic pulse can also be simulated in good agreement with the experimental data.
the value of the objective function increased only by  <dig> %, when the degrees of freedom increased by 21% , to  <dig> ) showing that elimination of these four elasticities in fact improved the quality of the fit.
cross validation of the obtained model
in order to validate the estimated elasticities, we have performed a cross validation study, for which we have generated an independent dataset, consisting of steady state perturbation data, and we have predicted the effect of that perturbation with our estimated parameters, and compared the results.
since we already perturbed vin, vgapd and vatpase and included the resulting data in the identification dataset , we had to select another enzyme in the pathway as the target for independent steady state perturbation.
among the remaining enzymes, we chose to inhibit the activity of vpfk by 50%, because this enzyme assumes many interactions and is known to be a complex enzyme, which would be a challenge to cross validate for our network.
the comparison between the model prediction and the experimental results is given in figure  <dig>  the level of extracellular and intracellular glucose did not change, whereas the levels of g6p and polysaccharides increased considerably and the levels of pep, fdp, ethanol and glycerol decreased.
the linlog model predicted this steady state perturbation with maximal deviation of only  <dig> % in g6p.
we conclude that the linlog model with the estimated parameters performed satisfactorily for the test case considered.
calculation of the systemic properties
after the estimation of the elasticities, we have calculated the flux control coefficients  using equation .
the comparison of these with the theoretical control coefficients of the main flux  is presented in figure  <dig>  the flux is controlled mainly by three enzymes: the hexose transporter , and to a lesser extent by the phosphofructokinase , and the atpase .
these findings are in qualitative agreement with the literature  <cit> .
it is noticeable that the elimination of the four elasticities did not change these results significantly; the same three reactions still have the main control of the glycolytic flux.
the concentration control coefficients  are also calculated using equation .
figure  <dig> gives the comparison of the cx0s of the two branch point metabolites, g6p and fdp.
the levels of both metabolites can be increased by increasing the activity of the hexose transporter.
increasing the activity of eatpase will decrease the level of g6p, whereas it increases the level of fdp.
increasing the activity of epfk will decrease and increase the levels of g6p and fdp respectively.
lastly, an increase in the activity of egapd will result in a slight increase in g6p level and a decrease in fdp level.
these results are also in agreement with the previous studies.
as expected, the model reduction did not have a significant effect on the estimated control coefficients.
discussion
large mathematical models are needed for finding the targets for engineering realistic metabolic systems.
currently, the available large models of several organisms e.g.
saccharomyces cerevisiae are mostly stoichiometric models.
although these models are highly relevant for knockout studies and determining the capabilities of the network considered, they have some limitations, e.g., they cannot predict the time courses, do not give insights in molecular events and are based on the "assumed goals" of the cells.
hence, for driving engineering interventions, kinetic models are needed.
here, we have presented a method to estimate the elasticities using the linlog kinetic format, from a combination of steady state and dynamic pulse experiments.
in the past, there have been efforts to extract mca parameters from transient metabolite data.
liao and delgado presented a method to obtain the control coefficients from such data  <cit> .
later, it was shown by elde and zacchi through mc simulations that this method is highly sensitive to noise  <cit> .
moreover a full dynamic kinetic model is not obtained.
here, we have described an estimation procedure which directly yields the elasticities as parameters instead of control coefficients, hence allowing the construction of a full kinetic model.
subsequently, the control coefficients can be obtained from the summation and connectivity theorems.
such a full kinetic model can also be used to simulate the effect of different perturbations on metabolism of a microorganism in a fermentor.
the method presented here assumes the availability of dynamic perturbation data.
although the list of metabolites is growing, not all of the metabolites in the cell can be measured.
in a recent study, wang and coworkers described an extension of the mca, under uncertainty  <cit>  in which the authors described a framework to calculate the control coefficients when either there are no measurements on metabolites, or the available measurements are subject to high uncertainty.
their proposed framework can be considered as complementary to the method presented in this paper.
in a case where accurate measurements are available, the present paper provides a mathematical approach to obtain the elasticities as kinetic parameters, from which one can calculate the control coefficients.
when compared to previous studies where linlog kinetics have been applied in combination with highly idealized linear pathways or small networks with branch points, the model in this work represents a more realistic case.
the model used in this study represents an intermediate size system with branch points and conserved moieties.
the proposed method can easily be extended to larger networks that are needed not only to understand the functioning of living cells, but also to infer engineering applications of e.g.
microorganisms.
in the example considered in this paper, we have assumed that the zero entries of the elasticity matrix are known.
this is a reasonable assumption, since numerous enzymes in the primary metabolism of many organisms are extensively studied and there is a dedicated public compendium for such information  <cit> , allowing the use of a priori knowledge on elasticities which are zero.
in the absence of any a priori knowledge on the zero elasticities, the elasticities have to be estimated from the full elasticity matrix, for which vastly more experimental data are needed.
in the mca literature, there are several attempts to determine the different perturbations needed in order to determine all elasticities  <cit> .
however, one needs to consider that an enzyme will never be affected by all metabolites; there is a physical limit for e.g.
the maximum number of binding sites of each enzyme in a biological network.
from the results presented here, it has been shown that, despite the rich information content of the data obtained from dynamic experiments, not every elasticity of the network could be correctly estimated.
this problem of parameter identifiability would be much more pronounced in a case where all  interactions are taken into account, caused by the combinatorial explosion of number of parameters to be estimated.
in the current work, in order to resolve some of the interactions which could not be resolved from the dynamic data, steady state enzyme perturbations have been introduced.
in order to get full kinetic models of microorganisms, accurate measurements of metabolite concentrations resulting from independent perturbations are needed, such as presented in  <cit>  and  <cit> .
at this point it is important to state that the property of the linlog kinetic format, that the rate equations are linear in the kinetic parameters, allows simultaneous use of alternative datasets by concatenating them in one parameter estimation scheme, i.e.
it is straightforward to extend the data matrix y in equation  and the vector χsim in equation , with the data from alternative dynamic and steady state perturbation experiments.
in addition, the decoupling of the parameters is also immediate i.e.
we can isolate poorly identifiable elasticities and estimate the remaining ones accurately and thereby reduce the problem.
it is worth mentioning that the obtained parameters  are true kinetic parameters that reflect the properties of the enzymes with respect to the corresponding metabolites.
they can be assumed to remain invariant as long as the enzyme keeps its properties in response to the changes in the environmental conditions.
additionally, since the linlog kinetics provides an approximation to the actual rate of the corresponding reaction over certain interval, the elasticity parameter that performs best may be different than the theoretical value of the corresponding elasticity.
this should be kept in mind when comparing some of the less well determined elasticity values to the theoretical ones.
CONCLUSIONS
constructing dynamic models of metabolic reaction networks under in vivo conditions using data obtained from perturbation experiments remains still a challenging problem in the area of systems biology.
in this contribution, we presented a method which allows the determination of the elasticities directly from experimental data from typical dynamic and/or steady state perturbation experiments.
these elasticities allow the reconstruction of the full kinetic model of the glycolysis of saccharomyces cerevisiae, and the determination of the control coefficients.
we further show by a posteriori parameter identifiability analysis that a subset of elasticities could not be identified using dynamic perturbation data only.
introduction of additional experimental information, i.e.
steady state experiments, solved this parameter identification problem.
considering the description of the dynamics, results of the cross validation studies and the final values of the elasticities, we conclude that linlog kinetics, although being an approximate kinetic format, performs very satisfactorily for estimating elasticities from data obtained from dynamic simulation of a mechanistic model, that was realistic in terms of biochemical complexity , noise added and sampling frequency.

 metabolic activities of living cells are accomplished by a well regulated, highly coupled network of numerous enzyme catalyzed reactions and selective membrane transport systems. to engineer such systems, enzymatic, transport and regulatory functions of the cells are manipulated via the use of recombinant dna technology. within the purpose of metabolic engineering, i.e. rational redesign of the metabolic systems, the highly relevant question of which perturbation should be applied in order to increase the productivity of the microorganism is addressed. the answer to this question requires information on both the regulatory level and the metabolic level. ter kuile and westerhoff showed by their hierarchical control analysis that the pathway flux is rarely controlled solely by gene expression, but that metabolite levels are also relevant. in this context, we focus on the metabolic reaction network level. modeling metabolic reaction systems is usually based on stoichiometric sometimes followed by kinetic modeling. one of the initial steps in the modeling of metabolic reaction networks is to determine the structure and steady state characteristics of a given network using stoichiometric information alone. steady state models describe time invariant fluxes, gathered from steady state experiments; hence they reflect the structural characteristics of the system. metabolic flux analysis and metabolic network analysis were developed as powerful tools to analyze such flux data. at steady state, the mass balances over the metabolites in the metabolic network yield a set of linear relations between the metabolic fluxes which can be expressed as: s·v =     where s is the stoichiometric matrix and v is the vector of metabolic fluxes, where m is the number of balanced metabolites and r is the number of fluxes. here, the system is in most of the cases underdetermined so that there are an infinite number of possible solutions. the realized solution depends on the kinetic properties of the involved reactions; this information is seldom known. to bypass the need of information on kinetics of individual reactions, alternative mathematical approaches have been proposed in the past to obtain a unique solution. an example is the constraint based optimization approach which is based on assumed optimality criteria, e.g. maximum growth, given biochemical, thermodynamic and irreversibility constraints and maximal reaction rates. later, segre et al. proposed the optimality constraint that requires maximization of biomass formation while minimization of metabolic adjustment in order to obtain a unique flux distribution of a mutant strain. in their approach, they defend their optimality criterion that a knock-out mutant strain would optimize its biomass production rate by changing minimally its metabolic fluxes from the wild type strain. an alternative metabolic modeling framework which uses a fitness function is the cybernetic approach. this approach assumes that an organism is an optimal strategist in utilizing all available sources with maximum efficiency. the expression and activity of the enzymes that catalyze network functionality are regulated by cybernetic control variables obtained from the solution of a constrained optimization problem. despite a number of successful applications especially in mixed substrate and prediction of knock-out lethality, all stoichiometric modeling approaches have their limitations, e.g. they can not predict time courses of the cellular processes, are based on "assumed goals" of the cell, and do not give insight in molecular events occurring in the cells, since information about the kinetic properties of the individual enzymes are not required. in the light of the above arguments, it is apparent that, to advance our understanding of the functioning of living cells a systems biology approach is needed, whereby the use of dynamic mathematical models of metabolic reaction networks to describe the complex kinetic behavior and interactions is becoming increasingly relevant. the kinetic behavior of many important enzymes occurring in metabolic networks have been studied extensively, however these studies have generally been performed under non-physiological conditions in test tubes, and therefore the applicability of these results to the in vivo metabolism is doubtful. teusink et al. showed that discrepancies exist between the in vivo measured changes of the concentrations of the glycolytic metabolites and their estimates using models based on mechanistic rate equations and in vitro determined parameters. this basic problem invalidates detailed models of metabolism containing kinetic parameters which have been determined in vitro. therefore, it is preferred to base the kinetic analysis of metabolic networks on in vivo studies of intact cellular networks. these in vivo studies are based on steady state and/or dynamic perturbations of a metabolic network starting from a reference steady state that is defined by its fluxes, enzyme activities, metabolite levels, extracellular concentrations, and kinetic parameters. it is also important to notice that the number of parameters that is typically involved in the traditional mechanistic equations is very large which causes identifiability problems due to parameter insensitivity. despite the information richness of data obtained from dynamic perturbations in in vivo experiments, there is a limit for identification of the parameters. for these reasons, it seems justified not to limit ourselves to the available complex mechanistic enzyme kinetics which must also be considered as approximations of the true in vivo behavior. approximative enzyme kinetic formats, which contain much fewer kinetic parameters are therefore of general interest. an overview of different approximative kinetic formats used in metabolic network modeling is given by heijnen. one of the proposed formats is, linlog kinetics, which has been introduced for modeling of in vivo kinetics and for metabolic redesign, and shown to have a good approximation quality, standardized format and relatively few parameters. in linlog kinetics, all the rate equations have the same mathematical structure in which the relation between rates and enzyme levels is proportional, while for metabolite levels, a linear sum of logarithmic concentration terms is proposed. all variables are considered relative to a reference steady state ). the linlog approximation is valid in the neighborhood of the reference state ), but quite large changes of metabolite concentrations, enzyme levels and fluxes are allowed. the parameters in the kinetic equations are the same scaled elasticities mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadaqadaqaaggaciab=v7alnaadaaaleaacqwg4baeaeaacqwg2bgdaagccqgh9aqpdawcaaqaaiabdiha4naacaaaleqabagaegimaadaaagcbagaemosao0aawbaasqabeaacqaiwaamaaaaaowaaqgaaeaadawcaaqaaiabgkgi2kabdaha2bqaaiabgkgi2kabdiha4baaaiaawia7amaabaaaleaacqwg4baecqgh9aqpcqwg4baedaahaaadbeqaaiabicdawaaaasqabaaakiaawicacaglpaaaaaa@45b3@ that are used in metabolic control analysis. it is important to note that the elasticity parameters appear in the model in a linear fashion. vj0=ee0+εcvln⁡)      mathtype@mtef@5@5@+=feaafiart1ev1aaatcvaufkttlearuwrp9mdh5mbpbiqv92aaexatlxbi9gbaebbnrfifhhdyfgasaach8aky=wiffydh8gipec8eeeu0xxdbba9frfj0=oqffea0dxdd9vqai=hguq8kuc9pgc9s8qqaq=dirpe0xb9q8qilsfr0=vr0=vr0dc8meaabaqaciaacagaaeqabaqabegadaaakeaadawcaaqaaiabdaha2bqaaiabdqeaknaacaaaleqabagaegimaadaaaaakiabg2da9maalaaabagaemyzaugabagaemyzau2aawbaasqabeaacqaiwaamaaaaaowaaewaaeaacqaixaqmcqghrawkiigacqwf1oqzdaqhaawcbagaemieaghabagaemodayhaaogagiibawmaeioba42aaewaaeaadawcaaqaaiabdiha4bqaaiabdiha4naacaaaleqabagaegimaadaaaaaaogaayjkaiaawmcaaiabgucariab=v7alnaadaaaleaacqwgjbwyaeaacqwg2bgdaagccyggsbabcqggubgbdaqadaqaamaalaaabagaem4yamgabagaem4yam2aawbaasqabeaacqaiwaamaaaaaagccagloagaayzkaaaacagloagaayzkaagaaczcaiaaxmaadaqadaqaaiabikdayagaayjkaiaawmcaaaaa@5758@ when the elasticities are known, a full dynamic model of the whole metabolic network can be set up using linlog kinetics. such a model allows in principle the calculation of control coefficients also under dynamic conditions. in linlog kinetics, the elasticities are the kinetic parameters represented in the elasticity matrix. from these, and a given network structure specified in the form of a stoichiometric matrix, the control coefficients for a reference condition can be calculated from the summation and connectivity relations developed in the framework of mca. also the change in control coefficients upon large changes in enzyme levels can be calculated. moreover, the linlog formulation enables the analytical solution of the mass balances for steady state metabolite and flux levels in the metabolic network, providing the solution of the metabolic redesign problem, i.e. determination of the optimal enzyme levels that maximize a certain flux while the total amount of enzyme and the metabolite levels are constrained. visser et al., reported a successful application of linlog kinetics in an in silico study that aims to solve this metabolic redesign problem. in order to determine the kinetic parameters of a model of a given in vivo metabolic system, in vivo perturbations of the complete metabolic network have to be performed. there are two main types of perturbations that can be imposed on the system: steady state and dynamic perturbations. in steady state perturbations, usually the enzyme activity of one of the reactions is changed by adding specific inhibitors or activators or by genetically changing the enzyme activity, resulting in a new steady state. in steady state perturbations an important problem has been addressed by kacser and burns, which is the determination of the set of reactions that has to be perturbed in steady state fashion, in order to be able to determine all elasticities for a given metabolic network which resulted in their "double modulation". they showed that for a simple linear chain of reaction, perturbation of the activities of the first and last enzymes allows determination of the elasticities of all enzymes of the chain under the condition that each enzyme is only responsive to its substrate and product, which rules out the fact that feedback loops are present. the theoretical basis is presented in later studies which showed that determination of the elasticities for any enzyme in such a simple chain requires two perturbations, one upstream and one downstream of the enzyme concerned. giersch and cornish-bowden extended the double perturbation approach to more complex pathways containing branch points, regulatory loops, and conserved moieties and they provided guidelines to list the possible reactions to be modulated in order to determine the elasticities for arbitrary metabolic networks. from the obtained list a minimum set of steady state perturbations, to be imposed on a specific network, can be chosen. as an alternative to the analysis of giersch and cornish-bowden, hofmeyr and cornish-bowden offered co-response analysis, to identify the mono-functional units that respond together to any perturbation applied. these mono-functional units have to be dissected in order to determine the elasticities belonging to these groups. linlog kinetics has been successfully applied to estimate the elasticity parameters for a linear pathway from sets of steady state metabolite concentrations and enzyme activities. additionally, heijnen et al. proposed a method to obtain flux control coefficients around a branch point, from large enzyme perturbation leading to large steady state flux perturbations, using the linlog kinetic format. this approach allows obtaining explicit solutions for steady state flux and metabolite levels as a function of large changes in enzyme activities. an alternative to steady state perturbations are dynamic perturbations, in which the system, being initially at the reference state, is disturbed to create time dependent data of transient metabolite levels. the dynamic metabolite profiles are typically obtained as a series of snapshots in time. a recent example of these are in vivo measurements of a number of metabolites in rather dense time sequences of a few seconds or minutes, using'rapid sampling' methods with subsequent high precision metabolite measurement techniques. transient data are rich in information and allow determination of the time hierarchy of the different elements of the metabolic network and the causal relationships between the network elements. it is possible to exploit these data to estimate the parameters of a traditional full kinetic model. such a model can subsequently be used to calculate the values of the elasticities at a given reference steady state. alternatively, when the linlog kinetic format is used the elasticities can be directly estimated from the transient data. this was demonstrated by kresnowati et al. who estimated the elasticities of a small example network from transient concentration data, assuming that the dynamic fluxes are unknown. an important difference between steady state and dynamic perturbations is that in the former both the steady state enzyme levels and the metabolite concentrations have to be measured, whereas in rapid dynamic perturbation experiments only xi is required because the enzyme activities can be considered constant within a sufficiently short time window following the perturbation. note that flux data are required in both methods. however, fluxes are not independent variables as they follow from the measured intra- and extracellular metabolite concentrations and the proper mass balances. in this work, we present a method to estimate kinetic parameters, using linlog kinetics, using metabolome data obtained from steady state and dynamic perturbation experiments. as such, this can be considered as an extension of the work by visser and heijnen, wu et al., and kresnowati et al., who have provided the theoretical framework and a small practical application of the linlog kinetic format. we first generalize the notation to make it applicable to networks of arbitrary size and complexity. additionally, we address further the question of parameter identification from experimental data. we monitor the propagation of error throughout the proposed parameter estimation procedure and we determine the subset of identifiable parameters. to illustrate the proof of principle, we apply the presented theory on in silico generated data, with realistic experimental settings, to be able to compare the obtained results with the "known truth". this data is generated using a previously published glycolytic pathway model of saccharomyces cerevisiae.  in silico pulse experiments in these experiments, the steady state was perturbed by increasing the extracellular glucose concentration. when the external glucose concentration is increased, the intracellular glucose concentration also increases as the sugar is transported through facilitated diffusion. from figure 1a it can be seen that the g6p concentration increases around 6% in the first five seconds and then starts to decrease, as the atp concentration drops due to the phosphorylation of the glucose in the hexokinase reaction. atp starts to recover after about seconds due to increased increase atp production downstream of glycolysis. this is in agreement with the literature on the glucose uptake dynamics following a sudden increase in the external sugar concentration. following this atp fluctuation, the pep concentration also decreases and then recovers. in the first five minutes, the concentrations of the external metabolites increase because of their increased production. after this initial period of min, in which the external glucose concentration is nearly constant, the concentration of external glucose slowly decreases due to continuous consumption and wash-out from the reactor. the cells follow this glucose drop by dropping the levels of intracellular metabolites. for the extracellular products the washout from the reactor is larger than their production and hence their concentration also drops. after about minutes, the original steady state is restored. such a pulse experiment therefore always delivers a highly dynamic dataset, followed by a pseudo steady state dataset. parameter estimation the dynamic data obtained during the first five minutes of the glucose pulse experiment were used to estimate the elasticities via the linlog parameter estimation procedure outlined in the methods section. in figure 1a, the simulation of the same dynamic perturbation using the estimated elasticities is also given. in figure the estimated elasticities are compared with the theoretically calculated elasticities, derived from mechanistic rate equations at the reference state. although the linlog simulation results do not differ much from the noise-free experimental data generated with the mechanistic model, some estimated elasticities are far from the expected theoretical values. the difference is in some cases in magnitude, for other cases there is a sign contradiction. specifically, out of elasticities are predicted in very good agreement with the theoretical ones, estimated elasticities agree in sign, but the magnitudes differ by more than 50%, and elasticities are estimated with a sign contradiction. notice that in the case of e the 200% deviation is already expected. this elasticity explains the effect of g6p on the polysaccharide formation rate and in the mechanistic model this effect is explained with a hill type kinetic equation with a hill coefficient of, hence this rate is very sensitive to the changes in g6p. to mimic this behavior, e is estimated to be much higher than its theoretical value. the poor identification of the elasticities has two possible reasons: either the experimental design is poor or there is a structural problem resulting that some interactions in the network can not be resolved from the available information. from an experimentalist point of view, it can be seen from figure 1a that the changes in the intracellular metabolites during the first seconds are between 25–30% which can be easily detected with the current measurement techniques. the extracellular glucose and polysaccharides also change to a detectable extent, but the changes in the ethanol are only 2–3% which is hard to detect. to check if the poor experimental design causes problems in the parameter identification, we have altered the experimental design to create larger changes in these concentrations. first the biomass concentration was increased to gdw l- also at the same moment that the fermentor was pulsed with the glucose, the inflow and outflow of the fermentor was stopped so the operation was effectively switched from though-flow mode to batch mode, and the extracellular metabolites were not washed out anymore. this allowed more rapid accumulation of the secreted products resulting in much larger changes. also, the change in the extracellular glucose concentration was more pronounced because there was no further addition of glucose after the pulse. these changes did not result in considerable changes in the intracellular metabolite profiles; they follow the slight change in the external glucose profile. the results of the new experimental design are depicted in figure which represents the new experimental data where only the significantly changed external metabolite profiles are depicted, together with the linlog simulation using the estimated elasticities. however, no significant improvement was achieved in the identification of the elasticities compared to the initial experimental design. hence it can be concluded that the problem is not due to poor experimental design. a closer look at the estimated elasticities reveals that the elasticities belonging to the part of the metabolic network, between vhk and vgapd are correctly estimated, whereas the estimation of the elasticities belonging to the lower part of glycolysis and the uptake reaction is poor. this is mainly due to the fact that the information content of the metabolite data is insufficient to resolve the complex interactions of metabolites and enzymes. specifically, in the vpk reaction rate, the feed forward activation of fdp and the mass action effect of pep are assumed. the perturbation of the system via an increase in the external glucose concentration results in an increase in the vpk flux. however, this does not allow a separate determination as to what extent the change in vpk is due to mass action effect of pep or due to feed forward effect of fdp. this results in identifiability problems for the elasticity parameters that describe the ethanol production. the same argument is valid for the elasticities belonging to the glycerol producing reaction, since in the in silico network the glycerol production is assumed to be proportional with the ethanol production; hence the same combined effect is seen in the glycerol production, resulting in an identifiability problem for the latter triplet of elasticities in order to analyze the information content of the data, the fisher information matrix is calculated as described in. fim is an indicator of the information content with respect to the parameters and due to the parameter linearity in linlog kinetics fim is independent of the values of the kinetic parameters. it can be calculated as yty, y being the design matrix appearing in equation. the singular values of the fim hold information on the number of linear dependencies between the columns of the data matrix. in our system, we checked the singular values of fim and concluded that not all the singular values are equally significant, which shows that there are some colinearities within the data matrix. this is also evident from the condition number, which equals the ratio of the highest to the lowest eigenvalues of the matrix and is ideally in a completely uncorrelated case. in this case the condition number was calculated to be × this fact supports the previous argument that using only one single set of dynamic data is not sufficient to resolve all the elasticity values. in order to obtain better estimates of the elasticities, steady state perturbation data are introduced. the steady state enzyme perturbation data are generated as described in the methods section: in order to resolve the feed forward effect of fdp on the vpk reaction, the enzyme amount of the vgapd reaction is perturbed in addition to the vin and vatpase. using the steady state perturbation data in table and the calculation procedure described in the methods section, we have estimated the elasticities for the in silico network under study. after the estimation of the elasticities using steady state data only, we have also simulated the dynamic perturbation of figure 1a using linlog kinetics, using the estimated elasticities. the results are presented in figure figure 4a shows the experimental data and the linlog simulation of the same pulse, and figure 4b compares the elasticities estimated using the steady state data with the theoretical elasticities. as can be seen from figure 4b, using the steady state data, the elasticities are estimated in good agreement with the theoretical elasticities. however, as can be observed from figure 4a, the linlog model deviates significantly from the experimental data, so the estimated elasticities have to be further refined. having a good initial estimate for the linlog elasticities from the steady state perturbation data only, we used in the second step both the dynamic and steady state data for the parameter estimation procedure. the nonlinear regression procedure allows combining both available steady state and dynamic experimental data, so that all parameters can be accurately estimated. after the non-linear fit to the data, we have obtained the final estimation of the elasticities. the simulation of the pulse experiment with these final estimates, and the comparison of the elasticities are given in figure 5a and figure 5b respectively, which shows that most of the elasticities are correctly estimated and that the experimental data are well described. analysis of parameter identifiability and model reduction in the previous sections, we have seen that we could not determine all the parameters from dynamic perturbation data only and we have introduced the steady state perturbation data. after determination of all of the elasticities, using both types of perturbation data, we analyzed back the identifiability of each parameter under noise. in order to carry the identifiability and error propagation analysis, mc simulation was used. as described in the methods section, 10% relative error was added to the noise free data represented in figure 1a, and the non-linear estimation procedure was implemented using the elasticities estimated from the steady state data, as the initial guesses. after repeating this scheme times, we obtained a distribution for each of the elasticities from which we calculated the relative error for each elasticity as the standard deviation per mean of the corresponding elasticity. the adjustable parameter λ in equation determines the relative importance of the different sources of data. in our case, these are dynamic and steady state data. notice that during the mc simulations, noise is added only to dynamic data. for steady state perturbations, the number of data points that can be obtained is theoretically infinite. this implies that, considering the central limit theorem, the steady state data can, in principle, be considered as noise free. however, the aim of the mc simulations was to elucidate the potentially unidentifiable elasticities using dynamic perturbation data, and therefore we chose λ to be equal to zero. the other extreme would suppress the effect of noise and wouldn't lead to detection of poorly identifiable elasticities. the values in between, are up to the choice of the modeler, depending on how many data points have been obtained from a dynamic pulse experiment, the standard deviations of the measurements, etc. we have implemented different values for the value of λ, but the outcomes of the mc simulation, i.e. which elasticities can hardly be identified, did not change qualitatively. from the distribution of each elasticity, the p-value, indicating the probability that the actual elasticity is zero and that the estimated value was caused by the random error only, is calculated via t-test and the results are given in table. the elasticities e e e and e have a relative standard deviation higher than 200%, so these elasticities are not identifiable under the 10% noise present. therefore we conclude that these elasticities are clearly candidates for model reduction by assigning them to be zero. a closer look reveals that e represents the feed-forward effect of fdp on vpk. e and e represent the effect of atp on vhk and vgol respectively and e is the feedback inhibition of g6p on vin. notice that the elasticities of the reactions of vin, vpk and vgol were already badly identified, when we estimated the elasticities using dynamic data only. using the arguments stated above, we set the values of these four elasticities to zero, implemented the non-linear regression step, with the new sparser elasticity matrix, and estimated the elasticities using steady state and dynamic data. the results are presented in figure 6a, where the experimental data and the linlog simulation using the new set of elasticities are given and in figure 6b where the comparison of the estimated elasticities with the theoretical ones is presented. it is clear that with the new set of elasticities the dynamic pulse can also be simulated in good agreement with the experimental data. the value of the objective function increased only by %, when the degrees of freedom increased by 21%, to ) showing that elimination of these four elasticities in fact improved the quality of the fit. cross validation of the obtained model in order to validate the estimated elasticities, we have performed a cross validation study, for which we have generated an independent dataset, consisting of steady state perturbation data, and we have predicted the effect of that perturbation with our estimated parameters, and compared the results. since we already perturbed vin, vgapd and vatpase and included the resulting data in the identification dataset, we had to select another enzyme in the pathway as the target for independent steady state perturbation. among the remaining enzymes, we chose to inhibit the activity of vpfk by 50%, because this enzyme assumes many interactions and is known to be a complex enzyme, which would be a challenge to cross validate for our network. the comparison between the model prediction and the experimental results is given in figure the level of extracellular and intracellular glucose did not change, whereas the levels of g6p and polysaccharides increased considerably and the levels of pep, fdp, ethanol and glycerol decreased. the linlog model predicted this steady state perturbation with maximal deviation of only % in g6p. we conclude that the linlog model with the estimated parameters performed satisfactorily for the test case considered. calculation of the systemic properties after the estimation of the elasticities, we have calculated the flux control coefficients using equation. the comparison of these with the theoretical control coefficients of the main flux is presented in figure the flux is controlled mainly by three enzymes: the hexose transporter, and to a lesser extent by the phosphofructokinase, and the atpase. these findings are in qualitative agreement with the literature. it is noticeable that the elimination of the four elasticities did not change these results significantly; the same three reactions still have the main control of the glycolytic flux. the concentration control coefficients are also calculated using equation. figure gives the comparison of the cx0s of the two branch point metabolites, g6p and fdp. the levels of both metabolites can be increased by increasing the activity of the hexose transporter. increasing the activity of eatpase will decrease the level of g6p, whereas it increases the level of fdp. increasing the activity of epfk will decrease and increase the levels of g6p and fdp respectively. lastly, an increase in the activity of egapd will result in a slight increase in g6p level and a decrease in fdp level. these results are also in agreement with the previous studies. as expected, the model reduction did not have a significant effect on the estimated control coefficients. discussion large mathematical models are needed for finding the targets for engineering realistic metabolic systems. currently, the available large models of several organisms e.g. saccharomyces cerevisiae are mostly stoichiometric models. although these models are highly relevant for knockout studies and determining the capabilities of the network considered, they have some limitations, e.g., they cannot predict the time courses, do not give insights in molecular events and are based on the "assumed goals" of the cells. hence, for driving engineering interventions, kinetic models are needed. here, we have presented a method to estimate the elasticities using the linlog kinetic format, from a combination of steady state and dynamic pulse experiments. in the past, there have been efforts to extract mca parameters from transient metabolite data. liao and delgado presented a method to obtain the control coefficients from such data. later, it was shown by elde and zacchi through mc simulations that this method is highly sensitive to noise. moreover a full dynamic kinetic model is not obtained. here, we have described an estimation procedure which directly yields the elasticities as parameters instead of control coefficients, hence allowing the construction of a full kinetic model. subsequently, the control coefficients can be obtained from the summation and connectivity theorems. such a full kinetic model can also be used to simulate the effect of different perturbations on metabolism of a microorganism in a fermentor. the method presented here assumes the availability of dynamic perturbation data. although the list of metabolites is growing, not all of the metabolites in the cell can be measured. in a recent study, wang and coworkers described an extension of the mca, under uncertainty in which the authors described a framework to calculate the control coefficients when either there are no measurements on metabolites, or the available measurements are subject to high uncertainty. their proposed framework can be considered as complementary to the method presented in this paper. in a case where accurate measurements are available, the present paper provides a mathematical approach to obtain the elasticities as kinetic parameters, from which one can calculate the control coefficients. when compared to previous studies where linlog kinetics have been applied in combination with highly idealized linear pathways or small networks with branch points, the model in this work represents a more realistic case. the model used in this study represents an intermediate size system with branch points and conserved moieties. the proposed method can easily be extended to larger networks that are needed not only to understand the functioning of living cells, but also to infer engineering applications of e.g. microorganisms. in the example considered in this paper, we have assumed that the zero entries of the elasticity matrix are known. this is a reasonable assumption, since numerous enzymes in the primary metabolism of many organisms are extensively studied and there is a dedicated public compendium for such information, allowing the use of a priori knowledge on elasticities which are zero. in the absence of any a priori knowledge on the zero elasticities, the elasticities have to be estimated from the full elasticity matrix, for which vastly more experimental data are needed. in the mca literature, there are several attempts to determine the different perturbations needed in order to determine all elasticities. however, one needs to consider that an enzyme will never be affected by all metabolites; there is a physical limit for e.g. the maximum number of binding sites of each enzyme in a biological network. from the results presented here, it has been shown that, despite the rich information content of the data obtained from dynamic experiments, not every elasticity of the network could be correctly estimated. this problem of parameter identifiability would be much more pronounced in a case where all interactions are taken into account, caused by the combinatorial explosion of number of parameters to be estimated. in the current work, in order to resolve some of the interactions which could not be resolved from the dynamic data, steady state enzyme perturbations have been introduced. in order to get full kinetic models of microorganisms, accurate measurements of metabolite concentrations resulting from independent perturbations are needed, such as presented in and. at this point it is important to state that the property of the linlog kinetic format, that the rate equations are linear in the kinetic parameters, allows simultaneous use of alternative datasets by concatenating them in one parameter estimation scheme, i.e. it is straightforward to extend the data matrix y in equation and the vector χsim in equation, with the data from alternative dynamic and steady state perturbation experiments. in addition, the decoupling of the parameters is also immediate i.e. we can isolate poorly identifiable elasticities and estimate the remaining ones accurately and thereby reduce the problem. it is worth mentioning that the obtained parameters are true kinetic parameters that reflect the properties of the enzymes with respect to the corresponding metabolites. they can be assumed to remain invariant as long as the enzyme keeps its properties in response to the changes in the environmental conditions. additionally, since the linlog kinetics provides an approximation to the actual rate of the corresponding reaction over certain interval, the elasticity parameter that performs best may be different than the theoretical value of the corresponding elasticity. this should be kept in mind when comparing some of the less well determined elasticity values to the theoretical ones.  constructing dynamic models of metabolic reaction networks under in vivo conditions using data obtained from perturbation experiments remains still a challenging problem in the area of systems biology. in this contribution, we presented a method which allows the determination of the elasticities directly from experimental data from typical dynamic and/or steady state perturbation experiments. these elasticities allow the reconstruction of the full kinetic model of the glycolysis of saccharomyces cerevisiae, and the determination of the control coefficients. we further show by a posteriori parameter identifiability analysis that a subset of elasticities could not be identified using dynamic perturbation data only. introduction of additional experimental information, i.e. steady state experiments, solved this parameter identification problem. considering the description of the dynamics, results of the cross validation studies and the final values of the elasticities, we conclude that linlog kinetics, although being an approximate kinetic format, performs very satisfactorily for estimating elasticities from data obtained from dynamic simulation of a mechanistic model, that was realistic in terms of biochemical complexity, noise added and sampling frequency. 
dynamic modeling of metabolic reaction networks under in vivo conditions is a crucial step in order to obtain a better understanding of the functioning of living cells. the parameter estimation within the linlog kinetic framework as presented here allows the determination of the elasticities directly from experimental data from typical dynamic and/or steady state experiments. in this paper we present a method for estimation of the kinetic model parameters, which are equal to the elasticities defined in metabolic control analysis, from metabolite data obtained from dynamic as well as steady state perturbations, using the linlog kinetic format. these elasticities allow the reconstruction of the full kinetic model of saccharomyces cerevisiae, and the determination of the control coefficients. additionally, we address the question of parameter identifiability from dynamic perturbation data in the presence of noise. the method is illustrated using metabolite data generated with a dynamic model of the glycolytic pathway of saccharomyces cerevisiae based on mechanistic rate equations. mc simulations revealed that certain elasticities are potentially unidentifiable from dynamic data only. so far dynamic metabolic models generally have been based on mechanistic rate equations which often contain so many parameters that their identifiability from experimental data forms a serious problem. addition of steady state perturbation of enzyme activities solved this problem. finally, identifiable subset of parameters is determined using information on the standard deviations of the estimated elasticities through monte carlo simulations.
 dynamic modeling of metabolic reaction networks under in vivo conditions is a crucial step in order to obtain a better understanding of the functioning of living cells. so far dynamic metabolic models generally have been based on mechanistic rate equations which often contain so many parameters that their identifiability from experimental data forms a serious problem. recently, approximative rate equations, based on the linear logarithmic format have been proposed as a suitable alternative with fewer parameters.  in this paper we present a method for estimation of the kinetic model parameters, which are equal to the elasticities defined in metabolic control analysis, from metabolite data obtained from dynamic as well as steady state perturbations, using the linlog kinetic format. additionally, we address the question of parameter identifiability from dynamic perturbation data in the presence of noise. the method is illustrated using metabolite data generated with a dynamic model of the glycolytic pathway of saccharomyces cerevisiae based on mechanistic rate equations. elasticities are estimated from the generated data, which define the complete linlog kinetic model of the glycolysis. the effect of data noise on the accuracy of the estimated elasticities is presented. finally, identifiable subset of parameters is determined using information on the standard deviations of the estimated elasticities through monte carlo simulations.  the parameter estimation within the linlog kinetic framework as presented here allows the determination of the elasticities directly from experimental data from typical dynamic and/or steady state experiments. these elasticities allow the reconstruction of the full kinetic model of saccharomyces cerevisiae, and the determination of the control coefficients. mc simulations revealed that certain elasticities are potentially unidentifiable from dynamic data only. addition of steady state perturbation of enzyme activities solved this problem. 