Prior Work and Scope of Contribution
Mathematical and computational modeling of ME/CFS remains remarkably sparse relative to the disease’s complexity and the size of the affected population. This section surveys the existing literature honestly, identifies what has and has not been attempted, and positions the models developed in this part within that landscape. Transparency about the novelty—and the limitations—of the present work is essential for scientific integrity.
1 Existing Models
Four bodies of work constitute the substantive prior art in ME/CFS modeling.
Phair’s IDO metabolic trap. The most developed mechanistic model in the ME/CFS literature is Phair, Davis, and Kashi’s mathematical analysis of tryptophan metabolism (Phair, Davis, and Kashi 2019). This model demonstrates that substrate inhibition of indoleamine 2,3-dioxygenase 1 (IDO1), combined with non-functional IDO2 variants, produces a bistable system: one stable equilibrium at physiological tryptophan concentrations and a second pathological equilibrium at elevated concentrations. The model is genuinely novel in showing that a single-pathway bistability mechanism can explain both the triggering (transition to the pathological state during immune activation) and the persistence (stability of that state after the trigger resolves) of ME/CFS. It generates specific, testable predictions—including that recovery from the pathological state requires weeks even if tryptophan input is halted entirely. However, the model is narrowly scoped: it addresses one enzyme pathway and does not model energy metabolism, immune cell dynamics, neuroendocrine function, or the multi-system interactions that characterize the disease.
Broderick’s cytokine network analysis. Broderick and colleagues applied network graph analysis to cytokine co-expression patterns in CFS patients versus controls (Broderick et al. 2010). By measuring 16 plasma cytokines in 40 female CFS patients and 59 matched controls, they constructed mutual-information networks and identified topological differences: CFS networks were more hub-like, with attenuated Th1 and Th17 responses and altered NK cell signaling patterns. This work is computational but not dynamical—it characterizes network structure (which cytokines co-vary) rather than network dynamics (how cytokine concentrations evolve over time). It does not use ODEs, does not model temporal trajectories, and does not generate predictions about disease progression or treatment response. It remains the most rigorous network-level characterization of ME/CFS immune dysfunction.
Morris and Maes’s neuro-immune framework. Morris and Maes proposed a conceptual neuro-immune model describing how peripheral immune activation, oxidative and nitrosative stress, and mitochondrial damage interact to produce ME/CFS symptoms (Morris and Maes 2013). The model identifies the pathway from infection through immune activation, oxidative and nitrosative stress, and mitochondrial damage to post-exertional malaise. While comprehensive in scope and biologically well-grounded, it is a verbal/diagrammatic model, not a mathematical one. It does not specify rate equations, parameter values, or quantitative predictions. It identifies the structure of the problem that a mathematical model should formalize.
Li et al.’s genome-scale metabolic modeling. The most recent computational contribution is Li and colleagues’ application of genome-wide precision metabolic modeling (GPMM) to ME/CFS and Long COVID muscle tissue (Li et al. 2025). Using transcriptomic data and constraint-based flux balance analysis, they identified downregulation of the alanine/aspartate and arginine/proline metabolism pathways and proposed L-ornithine and L-aspartate (LOLA) supplementation as a candidate intervention. This work uses a genuine computational framework (constraint-based metabolic modeling), but it is static: flux balance analysis computes steady-state metabolic fluxes under optimality assumptions, without modeling temporal dynamics, feedback loops, or multi-system interactions. It does not address immune function, neuroendocrine regulation, or the dynamic phenomena (PEM, disease onset, relapse–remission) that define ME/CFS clinically.
2 Adjacent Computational Work
Two additional lines of work are relevant but do not constitute mechanistic modeling in the sense pursued here.
Xiong, Oh, and colleagues developed BioMapAI, a deep neural network trained on a four-year longitudinal multi-omics dataset from 249 participants (Xiong et al. 2025). BioMapAI integrates gut metagenomics, plasma metabolomics, immune cell profiling, and clinical symptoms to classify ME/CFS and identify disease-specific biomarkers. This is a powerful data-driven approach that provides systems-level insights, but it is fundamentally a pattern-recognition tool: it identifies statistical associations across omics layers without specifying the causal mechanisms, rate equations, or dynamical systems that generate those patterns. It cannot, by design, predict the consequences of perturbations (treatments, exercise challenges) that lie outside the training data distribution.
Similarly, metabolomic profiling studies—including the foundational work by Naviaux (Naviaux et al. 2016) and Germain et al. (Germain et al. 2020) —characterize the metabolic state of ME/CFS but do not model the dynamics that produce it. These studies are indispensable as data sources for parameterizing dynamical models but are not themselves models in the mathematical sense.
3 What This Part Contributes—and What It Does Not
The models developed in Chapters Energy Metabolism Models through Predictive Applications and Clinical Translation attempt something not, to our knowledge, previously published in the ME/CFS literature: a multi-system ODE framework that couples energy metabolism, immune dynamics, neuroendocrine regulation, autonomic control, coagulation, mast cell activation, connective tissue biomechanics, GI motility/SIBO dynamics, and epigenetic state dynamics into an integrated dynamical system of 64 state variables. Specifically, the contributions are:
- Mechanistic PEM model: To our knowledge, no prior work formulates post-exertional malaise as a dynamical system with explicit rate equations for ATP depletion, ROS generation, secondary immune activation, and delayed symptom amplification. The model extends to metabolic flexibility (Randle cycle, GPR81 feedback), carnitine shuttle dynamics, and mitochondrial quality control (fission/fusion/mitophagy/biogenesis) (Chapter Energy Metabolism Models).
- Extended immune and vascular models: Mast cell degranulation dynamics with self-amplification thresholds, coagulation cascade with microclot-mediated oxygen delivery reduction (predicting multiplicative VO2 impairment), endothelial dysfunction via eNOS uncoupling, and daratumumab response kinetics (Chapter Immune System Models).
- Neuroendocrine integration: BH4 three-way cofactor competition (TPH/TH/NOS) reveals coordinated serotonin/dopamine/NO deficits from a single bottleneck; cerebrovascular autoregulation models multiplicative CBF impairment; central sensitization is modeled as bidirectional feedback rather than passive readout; small fiber neuropathy couples to autonomic dysfunction (Chapter Neuroendocrine and Autonomic Models).
- Multi-system coupling: The integrated 64-variable model (Chapter Integrated Multi-System Models) formalizes bidirectional interactions—energy\(<->\)immune, neuro\(<->\)immune, HPA\(<->\)immune, cardiovascular\(<->\)metabolic, gut\(<->\)brain\(<->\)immune, autonomic\(->\)motility\(->\)SIBO\(->\)immune\(+\)energy, mast cell\(<->\)energy\(<->\)autonomic, coagulation\(<->\)oxygenation\(<->\)mitochondria—as explicit coupling terms. The model includes GI motility/SIBO dynamics with autonomic and mast cell drivers, connective tissue/EDS biomechanical coupling to POTS, and epigenetic state dynamics (methylation + acetylation as slow variables creating hysteresis). Phair modeled one pathway; Broderick characterized one network; the present work attempts to link these subsystems dynamically.
- Bifurcation and attractor analysis: Beyond Phair’s single-pathway bistability, the coupled system exhibits three distinct disease attractors (metabolic-dominant, immune-dominant, severe/locked), identified through bifurcation analysis. Separatrix topology determines which recovery paths are accessible from each attractor (Chapter Integrated Multi-System Models).
- Temporal evolution and early warning signals: Critical slowing down analysis yields wearable-detectable early warning signals for PEM crashes 24–48 h in advance. Hysteresis from epigenetic consolidation defines a 3–12 month intervention window after disease onset. Hopf bifurcation predicts 2–6 week symptom oscillation periods usable as proximity-to-recovery biomarkers. Kramers escape rate theory quantifies spontaneous recovery probability as a function of disease duration (Chapter~Temporal Evolution and Disease Trajectories).
- Clinical translation via control and network theory: Pacing is formulated as optimal control (Pontryagin maximum principle, Hamilton–Jacobi–Bellman for stochastic extension), yielding state-dependent rules and quantitative safety margins. Global sensitivity analysis (Sobol indices) ranks drug targets across the population and reveals subtype-specific target inversion. Network controllability analysis predicts that 4–6 simultaneous interventions are structurally required for full system control—providing a mathematical explanation for monotherapy failure. A synergy matrix identifies specific treatment combinations and antagonisms. Virtual population simulation enables in silico trial design with enrichment optimization (Chapter~Predictive Applications and Clinical Translation).
These contributions must be understood in the context of severe limitations.
No computational implementation. The models are presented as systems of equations, not as validated simulations. No numerical solutions, parameter fitting to patient data, or sensitivity analyses are reported. The models are a theoretical framework, not computational results. Whether the integrated system actually produces the predicted emergent behaviors (bistability, limit cycles, PEM dynamics) remains to be demonstrated through numerical exploration.
Standard mathematical techniques. The individual subsystem models use well-established biochemical modeling methods—Michaelis–Menten kinetics, Hill functions, compartmental ODE systems. The novelty lies in the application to ME/CFS and the inter-system integration, not in the mathematical methods themselves.
Parameter uncertainty. Most model parameters are constrained by in vitro studies on isolated cells (primarily lymphocytes and muscle biopsies). Whether these values accurately represent in vivo function across tissues is uncertain. The 64-variable integrated model has on the order of 220–270 parameters; current ME/CFS datasets cannot robustly constrain this many parameters simultaneously.
Validation gap. No model in this part has been validated against independent patient data. The perturbation response analyses (Chapter Integrated Multi-System Models, Section Whole-Body Systems Model) describe expected qualitative behaviors but do not constitute quantitative validation. Bridging this gap requires the large-scale longitudinal datasets currently under development (e.g., DecodeME (DecodeME Consortium, Ponting, et al. 2025)).
Scope relative to biology. Even at 64 state variables, the model omits substantial biology: microRNA networks, detailed neurotransmitter receptor pharmacology, intracellular signaling cascades (e.g., NF-\(\kappa\)B, JAK-STAT), and tissue-specific heterogeneity. The model is a simplification designed for tractability, not a comprehensive digital twin.
In summary, the present work is the most comprehensive attempt at dynamical systems modeling of ME/CFS to date—extending from single-subsystem ODEs through a 64-variable integrated model to bifurcation analysis, optimal control theory, and network controllability—but it remains a theoretical framework awaiting computational implementation and empirical validation. It should be read as a formalization of hypotheses—translating the verbal models of Morris and Maes (Morris and Maes 2013) and others into the language of differential equations, and then extracting predictions that are only accessible through the mathematical formalism—rather than as established quantitative science.