A persistent obstacle for constructing kinetic models of metabolism is uncertainty in the kinetic properties of enzymes. Currently, available methods for building kinetic models can cope indirectly with uncertainties by integrating data from different biological levels and origins into models. In this study, we use the recently proposed computational approach iSCHRUNK (in Silico Approach to Characterization and Reduction of Uncertainty in the Kinetic Models), which combines Monte Carlo parameter sampling methods and machine learning techniques, in the context of Bayesian inference. Monte Carlo parameter sampling methods allow us to exploit synergies between different data sources and generate a population of kinetic models that are consistent with the available data and physicochemical laws. The machine learning allows us to data-mine the a priori generated kinetic parameters together with the integrated datasets and derive posterior distributions of kinetic parameters consistent with the observed physiology. In this work, we used iSCHRUNK to address a design question: can we identify which are the kinetic parameters and what are their values that give rise to a desired metabolic behavior? Such information is important for a wide variety of studies ranging from biotechnology to medicine. To illustrate the proposed methodology, we performed Metabolic Control Analysis, computed the flux control coefficients of the xylose uptake (XTR), and identified parameters that ensure a rate improvement of XTR in a glucose-xylose co-utilizing S. cerevisiae strain. Our results indicate that only three kinetic parameters need to be accurately characterized to describe the studied physiology, and ultimately to design and control the desired responses of the metabolism. This framework paves the way for a new generation of methods that will systematically integrate the wealth of available omics data and efficiently extract the information necessary for metabolic engineering and synthetic biology decisions.
Kinetic models are one of the cornerstones of rational metabolic engineering as they allow us to capture the dynamic behavior of metabolism and to predict dynamic responses of living organisms to genetic and environmental changes. With reliable kinetic models, metabolic engineering and synthetic biology strategies for improvement of yield, titer, and productivity of the desired biochemical can be devised and tested in silico [1]. The scientific community has acknowledged the utility and potential of kinetic models, and efforts towards building large- and genome-scale kinetic models were recently intensified [2–9]. Nevertheless, the development of these models is still facing challenges, such as partial experimental observations and large uncertainties in available data [10–12]. The major difficulty in determining parameters of kinetic models are uncertainties associated with: (i) flux values and directionalities [13–16]; (ii) metabolite concentration levels and thermodynamic properties [13–16]; and (iii) kinetic properties of enzymes [2,17]. As a result of interactions of metabolite concentrations and metabolic fluxes through thermodynamics and kinetics, these uncertainties make parameter estimation difficult. Quantifying these uncertainties and determining how they propagate to the parameter space is essential for identification of parameters that should be measured or estimated to reduce the uncertainty in the output quantities such as time evolution of metabolites or control coefficients [18,19]. In biological systems, large uncertainties and partial experimental data commonly result in a population instead of in a unique set of parameter values that could describe the experimental observations. Such population of parameter sets is typically computed using Monte Carlo sampling techniques [3–5,8, 9,11,20–28]. However, the problem is when certain properties differ among models in a model population. For example, one such property is flux control coefficients (FCCs) [18,19,29]. In [30], we used the ORACLE (Optimization and Risk Analysis of Complex Living Entities) framework [3,4, 8,10,11,31,32] to compute a population of kinetic models along with the corresponding flux control coefficients with the aim of improving xylose uptake rate (XTR) of a glucose-xylose co-utilizing S. cerevisiae strain. We have found that in the same population of models that are consistent with the observed physiology FCCs can be different due to lack of data about kinetic parameters. This can lead to erroneous or conflicting conclusions and decisions about the system in metabolic engineering and synthetic biology studies. In this contribution, to resolve such issues, we propose to formulate these problems as parameter classification: identify which of the parameters, if any, should be constrained so that the values of studied properties, such as FCCs, are in predefined ranges. For this purpose, we extended the capabilities of iSCHRUNK (in Silico Approach to Characterization and Reduction of Uncertainty in the Kinetic Models), a recently introduced machine learning approach that characterizes uncertainties in parameters of kinetic models, and identifies accurate and narrow ranges of parameters that can describe a studied physiological state [17]. In iSCHRUNK, machine learning is combined with methods that generate populations of kinetic models [3–5,8, 9,11,20–28] to data-mine the integrated data and observed physiology together with the kinetic parameters. The extended iSCHRUNK workflow is amenable for identifying parameters that give rise to a wide variety of properties of metabolic responses. The identified parameters can further be refined in an iterative way using the stratified sampling. Moreover, a set of improvements in the parameter classification procedure was introduced to improve the classification accuracy and to allow for dealing with uncertainties in alternative physiologies, e. g. , when multiple metabolite concentrations vectors are consistent with the observed physiology. As an illustration of the capabilities of the extended iSCHRUNK, we identified the enzymes and their kinetic parameters that determine consistent FCC values related to XTR. Our results showed that by constraining only three parameters, corresponding to xylose reductase (XRI) and ATP synthase (ASN), consistent FCCs can be obtained for models computed around multiple steady-state metabolite concentrations. We further showed how the parameter classification can be improved to more accurately identify the parameter subspaces that lead to well-determined model properties.
In [30], we analyzed the improvement of the xylose uptake rate (XTR) during mixed glucose-xylose utilization in a recombinant Saccharomyces cerevisiae strain. Here, we revisited that study and built the kinetic model of S. cerevisiae metabolic network around the reference steady-state of metabolic fluxes and metabolite concentrations (Methods). The model contains 258 parameters and describes 102 reactions and 96 intracellular metabolites distributed over cytosol, mitochondria and extracellular environment. The experimentally determined values of kinetic parameters were missing, and the analyzed system was underdetermined, i. e. , we had 102+96 computed values for steady-state fluxes and metabolite concentrations versus 258 unknown parameters. This meant that a multitude of parameter sets could reproduce the observed physiology, and we used the ORACLE framework that employs Monte Carlo sampling to generate a population of 200’000 kinetic models. We computed the flux control coefficients (FCCs) of the metabolic network and used them to rank enzymes according to their control over XTR, i. e. , the highest ranked enzymes were the ones with the largest magnitude FCCs with respect to XTR. Among the top ranked enzymes, hexokinase (HXK), non-growth associated maintenance (ATPM), and NADPH reductase (NDR) had ambiguous control over XTR (Fig 1A). The distributions of the control coefficients of XTR with respect to HXK, ATPM and NDR (CHXKXTR, CATPMXTR, and CNDRXTR, respectively) were extensively spread around zero, and we could not deduce with certainty whether the control of these enzymes over XTR was positive or negative. The population of control coefficients CHXKXTR was nearly symmetric around zero with a mean of 0. 005 and 47% of samples had negative values (Fig 1B). We split the population of kinetic models based on the sign of CHXKXTR, and we analyzed the two populations with a negative (Fig 1C, left) and a positive (Fig 1C, right) control of HXK over XTR. The split in the population did not have a substantial effect on the majority of the control coefficients. Interestingly, the exceptions were precisely the other enzymes with the ambiguous control over XTR, i. e. , ATPM and NDR, which exhibited a negative correlation with HXK (Fig 1C). This suggested that there were two distinct populations of kinetic models. The fact that models within these two populations have several common metabolic responses further implied that each of these two populations of models had distinct values of some kinetic parameters that determined such metabolic responses. We used the Classification and Regression Trees (CART) algorithm [33,34] to identify significant parameters that determine responses of XTR to changes in HXK activity. The CART algorithm partitions the parameter space into hyper-rectangles determined by the ranges of parameters that satisfy the studied property. Here, we used as parameters the degree of saturation of the enzyme active site, σA [10], because this quantity is constrained in a well-defined range between 0 and 1 (Methods), and the desired property was the negative control of HXK over XTR. The inputs of parameter classification were: (i) the information for each out of 200’000 parameter sets whether or not it gave rise to the negative control of HXK over XTR; and (ii) parameter values of 200’000 parameter sets. Subsequently, we will refer to hyper-rectangles computed by the CART algorithm as rules. To measure the performance of parameter classification we defined the performance index (PI), which quantifies a portion of parameter sets giving rise to the studied property. In this work, out of all parameter sets that satisfy rules (or a rule) inferred by parameter classification, PI quantifies how many of them are giving rise to the negative control of HXK over XTR. For example, within a population of models satisfying a rule, if 40% of models give rise to the negative control of HXK over XTR, then PI of this rule is 0. 4. To combine the distributions of top 3 parameters that ensure a high PI in an unbiased way, we performed another parameter classification (Methods). The parameter classification algorithm inferred 66 rules on these three parameters, and the top rule enclosed 9389 samples with PI of 0. 73 (S1 Table). The PI value of 0. 73 was close to the maximal PI value of 0. 78, which was computed for the rules formed with all parameters. As expected, the ranges of the three parameters defined by the top rule (Fig 4C) were consistent with the analysis presented in the previous section. We proceeded with the validation of the ranges of the top 3 parameters on a new population of models. We imposed the ranges of the top 3 parameters derived from the top rule of the parameter classification and generated a population of 100’000 models (Methods). We then computed the control coefficients of the top enzymes over XTR (Fig 4). The control coefficient CHXKXTR was distinctively negative with a mean value of -0. 09 (Fig 4C), and its distribution was clearly shifted toward negative values compared to that of the original population of models (Fig 1B). More than 72% of models had negative values of CHXKXTR compared to 47% in the original population of models. The value of PI of 0. 72 obtained from the validation set was strikingly close to the predicted value of 0. 73 from the second tree training. In Metabolic Control Analysis (MCA), it is considered that the control coefficients depend only on elasticities, however this holds only when the reactions are irreversible and there are no conserved moieties. It has been shown that metabolite concentrations affect displacements of reactions from thermodynamic equilibrium, which in turn influence the control over fluxes and concentrations in the network [3,16,32,37]. Therefore, when there is uncertainty in physiology, e. g. , when several alternative concentration profiles correspond to experimental observations, the control coefficients derived from the kinetic models computed for these concentration profiles can be significantly different. iSCHRUNK can resolve this kind of problems by identifying the parameter values that give rise to well-determined control coefficients of the metabolic network for multiple alternative physiologies. As an illustration, we analyzed three alternative physiologies characterized with three alternative concentration profiles (Reference, Extreme1 and Extreme2) and a common flux profile (Methods). We have undertaken to identify significant parameters that ensure a well-determined control over XTR for these physiologies. For this purpose, we constructed two populations of 200’000 kinetic models for the Extreme1 and Extreme2 physiology (Methods). Overall, together with 200’000 parameter sets computed previously for the reference physiology (Reference), we had 600’000 parameter sets for parameter classification. In the three populations of models, 47% (Reference), 46% (Extreme1) and 39% (Extreme2) of the models had a negative CHXKXTR. We found no rule with PI equal to 1 in performed parameter classification studies. This suggested that the parameter subspaces leading to a negative and a positive control of HXK over XTR were not distinctly separated. To improve the parameter classification for the problems where the separation between the classes is fuzzy, we propose to employ the k-nearest neighbors (k-NN) algorithm (Methods). The k-NN algorithm allows us to identify the parameter sets from one class that are surrounded by the parameter sets of the other class and reassign them to the latter. In the context of finding parameter values that give rise to a certain property, this means that the parameter classification algorithm will find only those parameter sets that are surrounded by a majority of the parameter sets of the same class. This way, the separation between the classes will be increased at the expense of neglecting parameter sets from the regions with a heavy overlap of the classes. We reconsidered the classification for parameters determining a negative control of HXK over XTR in the Reference case, and we applied the k-NN algorithm with k = 10 over the set of initial 200’000 parameters in order to find the surrounding for each of parameter sets, and to perform the reassignment (Methods). If in the group of 10 closest neighbors of a parameter set the percentage of parameter sets from the same class was less than a reassignment threshold, r, then that parameter set was reassigned. We performed two parameter classification studies for two different reassignment thresholds, r, of 30% and 50% (Methods). We found that as the reassignment threshold was increasing the tree training algorithm was inferring a smaller number of rules (73 for r = 30% versus 31 for r = 50%). Furthermore, the inferred rules were enclosing a smaller number of parameter sets for higher values of r, i. e. , for r = 30% and 50%, the top rules enclosed respectively 13427 and 1339 parameter sets (Fig 5 and S1 Table). In contrast, the obtained PIs, were higher for r = 50% than for r = 30% (Fig 5). For example, PI of the top rule for r = 50% was 0. 83, whereas the one for r = 30% was 0. 73 (Fig 5 and S1 Table). A comparison between the original method with preselection, which is identical to the reassignment method with r = 0% (corresponding to no reassignment), and the reassignment methods for r = 30% and 50% showed a general tendency of the latter for obtaining rules with improved PI and that enclose a smaller number of parameter sets (Fig 5). We also tested the reassignment procedure for parameters determining a positive control of HXK over XTR in the case of the reference metabolite concentration with k = 10 and r = 60%. The classification algorithm inferred 19 rules with PIs ranging from 0. 75 to 0. 90. The rules were defined by only 28 parameters (S1 Table). The top rule enclosed 1711 parameter sets with PI of 0. 90, and it was defined by 6 parameters. To validate the proposed improvement to the parameter classification, we imposed the distributions of the parameters defined by the top rules for the negative control case with r = 50%, and for the positive control case with r = 60% (S1 Table). We generated for each study a population of 100’000 models, and we computed the control coefficients in the network. In the case of negative control, the distribution of the control coefficient CHXKXTR was biased toward negative values with mean -0. 13 (Fig 6A and 6B). More than 79% of the computed control coefficients CHXKXTR were negative. (Table 4). Similarly, in the case of positive control, the distribution of the control coefficient CHXKXTR was shifted toward positive values with the mean of 0. 21 and a remarkable PI of 0. 89 (Fig 6D and 6F, and Table 4). For the negative and positive cases, the top rules were defined by 6 parameters each, where three parameters, σnadh_cXRI, σpi_mASN, and σt3p_cTPI, were common for both cases (Fig 6C and 6E). These three parameters were also ranked as the top 3 parameters in the parameter classification with the original algorithm (Fig 4). Moreover, the range of σnadh_cXRI was constrained toward low values for the negative control and toward high values for the positive control. In contrast, the parameters σpi_mASN and σt3p_cTPI were constrained toward high values for the negative control, and toward low values for the positive control. These patterns suggest that these three parameters are crucial for determining the sign of the control of HXK over XTR, whereas the remaining parameters, σdhap_cGPD1, σatp_mASN, and σxlt_cXRI for the negative control case, and σnad_cXRI, σatp_cATPM, and σnadp_cNDR for the positive control case, are likely having a minor effect on the PI. This result clearly demonstrated that the reassignment procedure allows for more precise identification of the subspaces leading to a desired control of HXK over XTR. We observed improvement of both PI and the mean CHXKXTR value compared to the results obtained with the unaltered parameter classification algorithm.
Machine learning methods [38–42] have found applications in a large number of biological and biomedical areas such as cancer research [43–45], population genetics [46,47], protein structure and function prediction and phylogenomic mapping [48–52], protein-protein interactions [53–55], medical imaging [56–60], gene expression and microarray data analysis [61–64], regulatory interactions [65,66], metabolic pathway dynamics [67], biomarker discovery and analysis of metabolomics and proteomics data [68–71]. However, the potential of these methods for detecting patterns in parameters of kinetic models of metabolism and uncovering hidden relationships between kinetic parameters, omics data, and observed phenotypes remained largely unexploited. Machine learning methods require large sets of training data for their successful application and methods for generating kinetic metabolic models that use Monte Carlo sampling offer an unprecedented opportunity for employing machine learning to advance our understanding of metabolic processes in cellular organisms. Kinetic models are usually built around a metabolic steady-state, which is characterized by the metabolite concentrations and metabolic fluxes, and the generated populations of kinetic parameters together with the observed steady-state data contain implicit information about the studied physiology. This information, if extracted systematically, can be used as guidance for the design of metabolic engineering and synthetic biology strategies that ensure the desired metabolic responses of studied organisms. In this work, we have extended iSCHRUNK functionalities to data-mine this information and systematically reduce uncertainties in the values of kinetic parameters that give rise to the desired metabolic behavior. As a demonstration, we reduced the uncertainties in the kinetic parameters that ensure that values of flux control coefficients remain within a pre-specified range. iSCHRUNK lends itself to a broad scope of applications ranging from sustainable production of biochemicals to medicine and regarding both the analysis and design of metabolism. It allows us to analyze the relationships between the inferred parameter ranges and the measurements acquired on the actual biological system, and, consequently, to create hypotheses regarding the operating states of enzymes and provide information about saturations of all enzymes in the network. This information is crucial for biotechnology studies where living cells need to be engineered for improved performance, or for drug discovery studies where, e. g. , we want to overproduce a compound that is toxic to a pathogen. The method can be applied not only to identify distributions of kinetic parameters but also to determine distributions of the metabolic fluxes and metabolite concentrations satisfying given requirements. It can also be used for guaranteeing both qualitative and quantitative features of metabolism, and several requirements can be combined simultaneously. For example, iSCHRUNK can be used to identify and quantify the parameters that maintain a redox potential while enforcing the desired level of yield and specific productivity of a compound of interest. Provided that the desired properties are biologically feasible, the method can be used to guarantee an arbitrary number of requirements. Finally, iSCHRUNK can be used to alleviate issues with high computational requirements of Monte Carlo sampling of kinetic parameters in large- and genome-scale metabolic networks. As the size of the models and complexity of studies increases, sampling a kinetic space becomes increasingly difficult and even intractable. However, iSCHRUNK allows us to identify relevant kinetic parameters that correspond to the observed physiology. The key finding of the current and previous studies [17] is that only a small set of parameters corresponding to a few enzymes is sufficient to characterize the observed physiology. Therefore, once we identify the most relevant parameters, it suffices to densely sample the identified parameters while fixing the remaining parameters at arbitrary feasible values. This way, iSCHRUNK dramatically reduces the sampling space, thus enabling computational analyses of large-scale and genome-scale dynamic metabolic systems.
The computational method for characterization and reduction of uncertainty, iSCHRUNK, was proposed in [17]. iSCHRUNK involves a set of successive computational procedures that can help us to ascertain and quantify the kinetic parameters that correspond to a given physiology. iSCHRUNK can be used with any method that generates populations of kinetic models describing given physiology such as ensemble modeling [24] or ORACLE [3,4, 8,10,11,31,32]. Here, we extended the original iSCHRUNK workflow (17) by an iterative loop that uses parameter classification to perform stratified sampling of the kinetic parameters, i. e. , it allows identifying refined sets of parameters that lead to the desired metabolic behavior (Fig 7). We used the extended iSCHRUNK to identify the distribution of kinetic parameters that determine the sign in ambiguous distributions of control coefficients as follows: We carried out the parameter classification in several steps (Fig 7B). We first removed from the consideration the parameters that were not affecting the control over the analyzed flux. We then used the CART algorithm with the preselected parameters for three populations of kinetic models where each population was computed with a different metabolite concentration vector (see Step II of the framework discussed above). In the third step, we ranked the parameters over three concentrations, and we chose the top parameters to continue. We next refined the distributions of the top parameters for each concentration individually, and we then used this information to determine the consistent distributions of top parameters over all concentrations. We detail the parameter classification steps below. In the cases when the space of parameter sets leading to a negative and the one leading to a positive control over analyzed quantities are overlapping, the separation between parameter classes is fuzzy. To enhance the separation between the classes, we propose here utilization of the k-nearest neighbors (k-NN) algorithm in the parameter classification as follows [87]. For each of the parameter vectors, we first assessed whether or not they were determining, e. g. , a negative control, and we assigned them to two distinct sets. The first set, SN, contained parameter vectors that gave rise to a negative control, whereas the second set, SP, contained the ones that gave rise to a non-negative control. We then ran the k-nearest neighbors (k-NN) algorithm, and for each parameter vector from the set SN, we computed how many out of its k-nearest neighbors belonged to the same set (SN). For each of these parameter vectors, if the percentage of k-nearest neighbors that belonged to the set SN was higher than a pre-specified reassignment threshold, r, we then retained that vector in the set SN. For instance, for r = 50%, if more than 50% of k-nearest neighbors of the analyzed parameter set belonged to the set SN, that parameter set remained in the set SN. Otherwise, we re-assigned that parameter vector to the set SP. With the proposed reassignment procedure, we emphasized the regions of the parameter space that have a higher proportion of parameter vectors belonging to the set SN. The reassignment procedure introduced two new parameters: the reassignment threshold, r, and the number of nearest neighbors, k. The values of r were chosen on the basis of the initial, unbiased, sampling that was performed in Step III. Specifically, from the initial sampling we could assess the average percentage of SN parameter vectors in the set of all vectors. We then set r to be a larger than the average percentage so that the parameter classification algorithm could identify the regions in the parameter space with the above than average proportion of SN vectors. Assuming that the parameter space was sampled uniformly, we use the parameter k to choose the larger or smaller part of the parameter space around the analyzed parameter vector for a possible reassignment. Very large values of k are not recommended as the reassignment procedure would consider the overall parameter space and no samples would be retained in the set SN as r is chosen to be larger than the average percent of SN vectors in the overall set of parameter vectors. Bayesian inference relies on use of Bayes theorem to compute the conditional distribution of a parameter vector θ given observed data x: p (θ|x) =p (x|θ) p (θ) p (x) where p (θ|x) is the posterior distribution of the parameters θ, p (θ) is the prior distribution of parameters, p (x|θ) is the likelihood, and p (x) is the evidence. In computing the posterior distribution p (θ|x), the evidence can be ignored as it represents a normalizing constant. It is often computationally prohibitive to explicitly evaluate the likelihood function and Approximate Bayesian Computation (ABC) methods are used for approximating this function by simulations [88]. For this type of studies, the ABC rejection algorithm [89] can be used as follows. First, the prior distribution of kinetic parameters is generated using the ORACLE framework or any other method that uses Monte Carlo sampling of uncertain parameters for constructing populations of kinetic models [3–5,8, 9,11,20–28]. The corresponding control coefficients are next computed, and the parameter classification algorithm is then used to discard parameter vectors from the prior that gave rise to ambiguous control over analyzed quantities. As a result, the retained samples are distributed according to the approximate posterior distribution of kinetic parameters that give rise to well-determined control over analyzed quantities. The simulations in this study were performed in MATLAB using an Apple MacPro Workstation with 2. 7 GHz 12-Core Intel Xeon E5 processor and 64 GB of RAM memory. The required time to generate a set 200’000 kinetic models was ~12. 5h, whereas one run of the parameter classification algorithm required several minutes.
Kinetic models are the most promising tool for understanding the complex dynamic behavior of living cells. The primary goal of kinetic models is to capture the properties of the metabolic networks as a whole, and thus we need large-scale models for dependable in silico analyses of metabolism. However, uncertainty in kinetic parameters impedes the development of kinetic models, and uncertainty levels increase with the model size. Tools that will address the issues with parameter uncertainty and that will be able to reduce the uncertainty propagation through the system are therefore needed. In this work, we applied a method called iSCHRUNK that combines parameter sampling and machine learning techniques to characterize the uncertainties and uncover intricate relationships between the parameters of kinetic models and the responses of the metabolic network. The proposed method allowed us to identify a small number of parameters that determine the responses in the network regardless of the values of other parameters. As a consequence, in future studies of metabolism, it will be sufficient to explore a reduced kinetic space, and more comprehensive analyses of large-scale and genome-scale metabolic networks will be computationally tractable.
medicine and health sciences
chemical compounds
applied mathematics
metabolic networks
enzymology
carbohydrates
organic compounds
simulation and modeling
xylose
algorithms
physiological parameters
metabolites
mathematics
network analysis
enzyme metabolism
pharmacology
enzyme kinetics
drug metabolism
enzyme chemistry
research and analysis methods
computer and information sciences
chemistry
pharmacokinetics
biochemistry
organic chemistry
physiology
monosaccharides
biology and life sciences
physical sciences
metabolism
Uncertainty reduction in biochemical kinetic models: Enforcing desired model properties