Kinetic Pathway Simulation for mGWAS Interpretation

Summary

Kodate, Sato, Hishinuma, et al. (2025) use a kinetic (differential-equation) model of the human liver folate cycle, with individual enzyme reaction rates perturbed to simulate the effect of a genetic variant, to interpret metabolome-genome-wide association study (mGWAS) results from the Tohoku Medical Megabank cohort. Systematically reducing each enzyme's reaction rate and comparing the resulting steady-state metabolite concentration changes to the direction of mGWAS association effects showed strong qualitative agreement (Kendall's τ=0.77, p=0.04) for variant-metabolite pairs with small mGWAS p-values, while also revealing metabolite fluctuations for gene-metabolite pairs that were not significant in the mGWAS — evidence these associations are real but underpowered rather than false negatives.

Method

1. Kinetic pathway model: the Nijhout et al. differential-equation model of the human liver cell folate cycle (sourced from BioModels), spanning cytosolic and mitochondrial compartments, with experimentally derived initial metabolite concentrations and enzyme reaction rates.
2. mGWAS: ~90 million SNVs tested against blood metabolite concentrations (NMR via BOLT-LMM; MS via GCTA linear regression) in a discovery and replication cohort from the Tohoku Medical Megabank Project (ToMMo); 35 genome-wide significant peaks were identified, 27 replicated, spanning 115 genes.
3. Reaction-rate perturbation: for genes in the folate pathway, the corresponding enzyme's reaction rate was reduced (e.g., by 70% for the well-characterized MTHFR A222V variant, rs1801133) to simulate the variant's known effect on enzyme activity, and the resulting steady-state metabolite concentration shifts were compared to the direction (not magnitude) of the mGWAS association betas.
4. Systematic enzyme screen: every enzyme in the model was independently perturbed in 10% reaction-rate steps down to 10% of baseline, generating a full map of enzyme-to-metabolite sensitivity, then cross-referenced against all mGWAS results in the pathway.

Findings

• MTHFR validation: reducing the modeled MTHFR reaction rate by 70% (matching the experimentally reported activity loss for rs1801133/MTHFR A222V) produced a 14.7% decrease in steady-state formate concentration, closely matching a previously reported 15% decrease in blood plasma formate for carriers of this variant.
• Directional concordance with mGWAS: for MTHFR and FTD (formyltetrahydrofolate dehydrogenase) variants associated with multiple metabolites (homocysteine, glycine, formate, serine, N,N-dimethylglycine) at suggestive significance, the direction of simulated concentration change matched the sign of the mGWAS effect size in nearly all cases; across all matched variant-metabolite pairs, Kendall's τ=0.77 (p=0.04).
• Recovering underpowered associations: the systematic enzyme screen found several gene-metabolite pairs with large simulated concentration shifts that were not significant in the mGWAS, suggesting real biological effects that current mGWAS sample sizes lack power to detect — a testable prediction for future, larger mGWAS.
• Enzyme classification by robustness to variation: enzymes were categorized by simulated impact breadth — some (like MTHFR and FTD) affect multiple metabolites and also carry multiple common non-synonymous variants in the population, while several other pathway enzymes had minimal simulated impact and correspondingly no or few common non-synonymous variants, consistent with these genes being under stronger evolutionary constraint against activity-altering variation.
• Theoretical grounding: the qualitative direction-only concordance (rather than precise quantitative matching) is explained via "buffering structure" theory from chemical reaction network theory — subnetworks that are structurally isolated from the rest of the network such that reaction-rate perturbations within them cannot propagate outside, meaning steady-state directional responses to a rate change are a robust, model-structure-derived prediction even when exact kinetic parameters are uncertain.

Significance

This offers a complementary, single-pathway mechanistic validation approach to whole-network genome-scale flux prediction (see Metabolic Flux Modulation of Genetic Risk, which instead predicts organ-wide stoichiometric flux from imputed transcript abundance across ~500,000 individuals). Where that framework trades pathway-level kinetic detail for genome-scale coverage, this approach trades genome-scale coverage for a detailed, experimentally-grounded kinetic model of one well-characterized pathway — the two approaches represent different points on the same detail-versus-coverage tradeoff in linking genetic variants to metabolic phenotypes.

See Also

• Metabolic Flux Modulation of Genetic Risk — a genome-scale, stoichiometric-flux alternative to this pathway-level kinetic simulation approach.
• Recon3D — a genome-scale (rather than single-pathway) human metabolic reconstruction.

Citations

[1] Kodate, S., Sato, M., Hishinuma, E., Kojima, K., Motoike, I.N., ToMMo Study Group, Koshiba, S., Yamamoto, M., Yamada, K.D., & Kinoshita, K. (2025). Simulating metabolic pathways to enhance interpretations of metabolome genome-wide association studies. Scientific Reports, 15, 1634.