Diverse metabolic model parameters generate similar methionine cycle dynamics

Abstract

Parameter estimation constitutes a major challenge in dynamic modeling of metabolic networks. Here we examine, via computational simulations, the influence of system nonlinearity and the nature of available data on the distribution and predictive capability of identified model parameters. Simulated methionine cycle metabolite concentration data (both with and without corresponding flux data) was inverted to identify model parameters consistent with it. Thousands of diverse parameter families were found to be consistent with the data to within moderate error, with most of the parameter values spanning over 1000-fold ranges irrespective of whether flux data was included. Due to strong correlations within the extracted parameter families, model predictions were generally reliable despite the broad ranges found for individual parameters. Inclusion of flux data, by strengthening these correlations, resulted in substantially more reliable flux predictions. These findings suggest that, despite the difficulty of extracting biochemically accurate model parameters from system level data, such data may nevertheless prove adequate for driving the development of predictive dynamic metabolic models.

Figure 1. Methionine cycle as modeled by Reed et al (2004)

Modeled metabolites are in bold within rectangles; enzymes are in ovals; metabolite influxes are in italics. Unmodeled metabolites are not shown. The conversion of methionine (Met) to S-adenosylmethione (AdoMet) is catalyzed by either methionine adenosyl transferase I (MATI) or methionine adenosyl transferase III (MATIII), which are isozymes derived from the same polypeptide chain. The former is feedback inhibited by AdoMet, while the latter is upregulated by it in the concentration ranges used in this study. AdoMet is demethylated by either glycine N-methyltransferase (GNMT) or a general AdoMet dependent methyltransferases (METH) to give S-adenosylhomocysteine (AdoHcy). Both enzymes are feedback inhibited by AdoHcy. AdoHcy is converted to homocysteine (Hcy) by AdoHcy hydrolase. Hcy can be methylated by Met synthase (MS) or by betaine:homocysteine methyltransferase (BHMT, inhibited by AdoHcy and AdoMet) to re-form Met. Alternatively, it can combine with serine via cystathionine β-synthase (CBS, up-regulated by AdoHcy and AdoMet) to exit the cycle.

Figure 2. Simulated methionine flux perturbation and dynamic metabolic responses

Data points used for parameter identification are surrounded by error bars of ±15% to simulate experimental noise. (A) Methionine influx (Metin). (B) Methionine concentration (Met). (C) S-Adenosylmethionine concentration (AdoMet). (D) S-Adenosylhomocysteine concentration (AdoHcy). (E) Homocysteine concentration (Hcy). (F) Flux from Met to AdoMet (V_MATI + V_MATIII). (G) Flux from AdoMet into AdoHcy (V_METH + V_GNMT). (H) Flux between AdoHcy and Hcy (V_AH). (I) Flux from Hcy to Met (V_MS+V_BHMT). (J) Flux from Hcy to cystathionine (V_CBS).

Figure 3. Histograms of parameter values that reproduce the dynamic data shown in Figure 2 within the error bars

The results are from inversion of the simulated concentration data only (blue histograms), of both concentration and flux data (green histograms), and of both concentration and flux data with additional biochemical constraints on the search ranges for selected parameters (red histograms; constraints are shown by black arrows surrounding the red peaks for K^MATI_m, K^MATI_i, K^GNMT_m, K^GNMT_i, and K^BHMT_m; see Table 2 for exact ranges). For all parameters, the X-axis spans from 10^-2 to 10⁵. Global search was conducted over this 10⁷-fold range, except for the indicated constraints mentioned above for the red histograms. The X-axis has been binned into 50 equally spaced segments of the parameter search range (on the log₁₀ scale). The Y-axis values represent the frequency with which the identified parameter values fell within each bin.

Figure 4. Correlations between parameter values

Symbols indicate parameter combinations that reproduce the concentration data shown in Figure 2. MATI (A) and MATIII (D) catalyze the same reaction, as do MS (B) and BHMT (E). The conversion catalyzed by GNMT (C) is also catalyzed by METH, but no correlation plot shown for METH, as it is characterized by only a single parameter value. (A, B, C) Enzymes with minor fluxes. The solid dark lines correspond to parameter combinations producing 10% of the overall reaction flux of the reference model. For these enzymes, the observed parameter values are restricted to combinations of V_max and K_m that ensure low flux (< 10% of total). (D, E) Flux-dominant enzymes, for which V_max and K_m are roughly linearly correlated (R² ≥ 0.95).

Figure 5. Predictions of models that contain parameter vectors identified by inversion of the simulated data shown in Figure 2

(A) Simulated methionine (Metin, black line) and homocysteine (Hcyin, blue line) influx perturbations used to generate the output shown in the other panels. (B) Metabolite concentration responses. (C, D) Homocysteine methylation reaction flux (V_MS + V_BHMT) response. In panels (B) – (D), blue lines show the output of the reference model; grey lines (N = 200) show the output of the models containing the newly identified parameter vectors. The parameter vectors used in (B) and (C) were identified by inversion of concentration data only, and in (D) by inversion of both concentration and flux data.