Bifurcation Analysis and Disease Subtypes
The integrated model’s nonlinear feedback structure implies the existence of multiple steady states (attractors) for certain parameter ranges. Bifurcation analysis—systematic exploration of how the number and character of attractors change as parameters vary—provides the most powerful application of the mathematical framework, yielding results that are inaccessible to verbal reasoning alone.
1 Attractors as an Emergent Consequence of the Biochemistry
An epistemological point warrants explicit discussion. The ODEs in Energy Metabolism Models through Integrated Multi-System Models were built from biochemical literature—enzyme kinetics, immune cell dynamics, mitochondrial biophysics—without targeting any particular dynamical behaviour. A modeller working from the same literature would make similar (though not necessarily identical) choices, because the key structural elements—saturating production kinetics, immune cell activation dependent on energy availability, bidirectional coupling—are standard biochemical building blocks, not exotic constructions selected for their dynamical properties.
This matters because the bistability derived below is a consequence of these biochemical choices, not a premise. The model was not reverse-engineered from a desired phase portrait. However, we must be precise about what “emergent” means here: it does not mean that every possible model of ME/CFS biochemistry will exhibit bistability. Specific modelling choices—particularly the Hill exponent in the immune activation term (see Remark below)—are necessary conditions for bistability. What we claim is weaker and more defensible: given the standard biochemical forms used throughout Energy Metabolism Models and Immune System Models, bistability follows without additional assumptions.
This section demonstrates that claim analytically, by deriving attractor existence from first principles on a reduced two-variable subsystem.
1.1 The Reduced Energy–Immune System
Consider only the two state variables that carry the dominant positive feedback: ATP concentration \(A \equiv [\text{ATP}]\) and activated immune cell count \(I \equiv N_a\). All other variables are treated as quasi-steady-state parameters. From energy immune coupling and immune energy feedback, the reduced system is:
\[ \begin{aligned} \frac{d A}{d t} &= P(A) - \beta_0 - J_\text{immune}(I) \\ \frac{d I}{d t} &= k_\text{act}^\text{eff}(A) \cdot I_0 - (k_\text{exh} + d_a) I \end{aligned} \tag{1}\]
where \(P(A)\) is the net ATP production rate—a unimodal (hump-shaped) function of \(A\). The shape follows directly from the PFK-1 kinetics in pfk1: at low ATP, AMP is high (via adenylate kinase equilibrium $ 2 + $), maximally activating glycolysis, but the production machinery itself requires minimal ATP (hexokinase priming, ion gradients), so net production rises with \(A\); at high ATP, product inhibition of PFK-1 and low AMP reduce glycolytic flux, while ATP synthase slows as ADP becomes scarce (atp synthase), so production falls. The result is \(P(A)\) rising from near zero, reaching a peak at intermediate \(A\), then declining—a consequence of the adenylate kinase feedback loop, not a modelling choice. The remaining terms: \(\beta_0\) is basal demand, \(J_\text{immune}(I) = e_a I\) is immune energy drain (linear in activated cells), \(k_\text{act}^\text{eff}(A) = k_{\text{act,0}} A^2 \\/ (K_\text{ATP}^2 + A^2)\) is the Hill-type effective activation rate (immune energy feedback), and \(I_0\) is the resting immune cell pool (treated as constant in the reduced model). The parameter \(e_a > 0\) is the per-cell energy cost of immune activation.
1.2 Nullcline Analysis
The fixed points of the system are the intersections of the two nullclines: the \(A\)-nullcline (\(d A\\/d t = 0\)) and the \(I\)-nullcline (\(d I\\/d t = 0\)).
\(I\)-nullcline. Setting \(d I\\/d t = 0\):
$ I^*(A) = $ {#eq-i-nullcline}
This is a sigmoidal (Hill) function of \(A\), increasing from 0 at \(A = 0\) to a plateau \(I_max = k_{\text{act,0}} I_0 \\/ (k_\text{exh} + d_a)\) as \(A -> \infty\). The sigmoidal shape is a consequence of the Hill kinetics governing ATP-dependent immune activation (immune energy feedback).
The cooperativity exponent \(n = 2\) in immune energy feedback is a modelling choice with consequences for bistability. With \(n = 1\) (standard Michaelis–Menten), the \(I\)-nullcline becomes a hyperbolic function lacking the inflection point that enables triple intersection with the \(A\)-nullcline; the system is then monostable for all parameter values. With \(n \geq 2\), the \(I\)-nullcline is sigmoidal and triple intersection becomes possible.
The choice of \(n = 2\) is motivated by immunometabolism: immune cell activation involves cooperative processes (receptor clustering, signalling threshold effects, metabolic reprogramming switches) that produce ultrasensitive responses (Ferrell 1996). T cell activation in particular exhibits a switch-like dependence on TCR signal strength, well-described by Hill kinetics with \(n = 1.5\)–$ 3$ Altan-Bonnet and Germain (2005) However, \(n = 1\) cannot be excluded a priori for all immune cell types, and the bistability result depends on this choice. The analysis below should therefore be read as: if ATP-dependent immune activation is cooperative (\(n \geq 2\)), then bistability follows from the feedback structure.
\(A\)-nullcline. Setting \(d A\\/d t = 0\):
$ I_^A (A) = $ {#eq-a-nullcline}
As established above, \(P(A)\) is unimodal: it rises at low \(A\) (where the production machinery activates as minimal ATP becomes available), peaks at intermediate \(A\), and declines at high \(A\) (where product inhibition of PFK-1 and ADP scarcity throttle production). For the ME/CFS regime, where Complex I is partially inhibited (\(\alpha_\text{CI} < 1\)), \(P(A)\) is vertically compressed (lower peak) and horizontally compressed (the peak shifts leftward). The \(A\)-nullcline \((P(A) - \beta_0)\\/e_a\) inherits this unimodal shape: it rises at low \(A\), peaks, and descends at high \(A\).
Intersection multiplicity condition. Bistability exists when the \(I\)-nullcline and the \(A\)-nullcline intersect three times. The \(I\)-nullcline (i nullcline) is strictly monotone increasing and sigmoidal. The \(A\)-nullcline (a nullcline) is unimodal (one local maximum, no other critical points). A strictly monotone curve and a unimodal curve can intersect at most three times: at most once on the ascending limb, at most once near the peak, and at most once on the descending limb. Three intersections occur when the following condition holds:
$ lr(|){A=A} > lr(|){A=A} $ {#eq-bistability-condition}
That is, at the inflection point of the \(A\)-nullcline, its slope exceeds the slope of the sigmoidal \(I\)-nullcline. Expanding: the \(I^*\) slope at inflection is \(k_{\text{act,0}} I_0 \\/ (2 (k_\text{exh} + d_a) K_\text{ATP})\); the \(A\)-nullcline slope depends on \(P'(A) \\/ e_a\). The condition reduces to:
$ > $ {#eq-bistability-inequality}
1.3 Parameter Regimes: Healthy vs.~ME/CFS
bistability inequality can be evaluated for each parameter regime using the values in Mathematical Model Details.
Healthy regime (\(\alpha_\text{CI} = 1.0\), \(k_\text{exh} = 0.05 \text{day}^{-1}\), \(k_\text{recov} = 0.10 \text{day}^{-1}\)). With full Complex I activity, \(P'(A)\) is steep and \(P(A) - \beta_0 > 0\) for a wide range of \(A\). However, healthy \(k_\text{exh}\) is low and \(k_\text{recov}\) is high, so the immune exhaustion term \((k_\text{exh} + d_a)\) in the denominator of \(I^*(A)\) is moderate. The \(I\)-nullcline plateau \(I_max\) is therefore limited: the immune system self-regulates before escalating to energy-depleting levels. Numerically, the \(A\)-nullcline lies above the \(I\)-nullcline for all physiologically plausible \(A\), and the two curves intersect once—yielding a single, stable healthy fixed point. The eigenvalues of the Jacobian at this point are:
$ _{1,2} = [(f_{A A} + g_{I I}) sqrt((f_{A A} - g_{I I})^2 + 4 f_{A I} g_{I A})] $ {#eq-healthy-eigenvalues}
where \(f_{A A} = \partial(dot(A))\\/\partial A < 0\) (ATP dynamics are self-stabilising at steady state), \(g_{I I} = -(k_\text{exh} + d_a) < 0\) (immune populations decay in the absence of activation), and the cross-terms \(f_{A I} = -e_a < 0\) and \(g_{I A} = k'_{\text{act,0}} I_0 > 0\). For healthy parameters, \(|f_{A A} g_{I I}| > |f_{A I} g_{I A}|\), so the discriminant is negative (complex eigenvalues) or both real parts are negative. The single fixed point is a stable node or stable spiral: the healthy state is an attractor.
ME/CFS regime (\(\alpha_\text{CI} = 0.65\), \(k_\text{exh} = 0.15 \text{day}^{-1}\), \(k_\text{recov} = 0.04 \text{day}^{-1}\)). With reduced Complex I activity, \(P(A)\) is compressed: \(P(A) - \beta_0\) is positive over a narrower range of \(A\) and the peak of \(P'(A)\) is lower—the \(A\)-nullcline drops. Simultaneously, the ME/CFS regime has elevated \(k_\text{exh}\). This increases the denominator \((k_\text{exh} + d_a)\) of \(I^*(A)\), which lowers the plateau \(I_max\)—on its own this would make bistability harder, not easier. However, the critical effect is that the depressed \(A\)-nullcline (from reduced \(\alpha_\text{CI}\)) now has a lower, narrower hump that intersects the (also lowered) \(I\)-nullcline in the three-intersection geometry. Whether bistability inequality is satisfied depends on the quantitative balance between the compression of \(P(A)\) and the lowering of \(I_max\); for the stated parameter values, the condition holds (verified numerically in Mathematical Model Details).
Under ME/CFS parameters (\(\alpha_\text{CI} = 0.65\), \(k_\text{exh} = 0.15 \text{day}^{-1}\), \(k_\text{recov} = 0.04 \text{day}^{-1}\)) and with cooperative immune activation (\(n = 2\), The Hill Exponent Is a Modelling Choice), the reduced system (reduced system) has exactly three fixed points: two stable fixed points (nodes or spirals) and one saddle.
Step 1: A positively invariant region exists. Since \(A\) and \(I\) represent physical concentrations, the positive quadrant \(A \geq 0\), \(I \geq 0\) is preserved by the dynamics (at \(A = 0\): \(dot(A) = P(0) - \beta_0 \geq 0\) when \(I = 0\); at \(I = 0\): \(dot(I) = k_\text{act}^\text{eff}(A) \cdot I_0 > 0\) for \(A > 0\)). For the upper bounds, choose \(A_max\) so that \(P(A_max) < \beta_0\) (guaranteed by PFK-1 product inhibition at high ATP) and let \(I_max = k_{\text{act,0}} I_0 \\/ (k_\text{exh} + d_a)\) (the \(I\)-nullcline plateau). Define the trapping region as the set bounded by \(I = 0\), \(A = A_max\), \(I = I_max\), and the \(A\)-nullcline curve \(I = (P(A) - \beta_0)\\/e_a\) on the left (which lies to the left of the ascending limb and satisfies \(dot(A) > 0\) in its interior). On the upper boundary \(I = I_max\): \(dot(I) < 0\) because immune clearance \((k_\text{exh} + d_a)I\) exceeds activation at the plateau. On \(A = A_max\): \(dot(A) = P(A_max) - \beta_0 - e_a I < 0\). On \(I = 0\): \(dot(I) > 0\). On the left boundary (the \(A\)-nullcline curve): by construction \(dot(A) = 0\) on this curve and \(dot(A) > 0\) immediately to its right, so trajectories cannot exit leftward. This trapping region \(\mathcal{R}\) is positively invariant and contains all three fixed points.
Step 2: Exactly three fixed points. Fixed points are intersections of the \(I\)-nullcline and \(A\)-nullcline. The \(I\)-nullcline (i nullcline) is strictly monotone increasing (derivative always positive for \(A > 0\)). The \(A\)-nullcline (a nullcline) is unimodal: it has exactly one local maximum (inherited from the unimodal \(P(A)\)) and is therefore composed of two monotone pieces—an ascending limb and a descending limb. On any interval where both curves are monotone in the same direction, they can intersect at most once (by the monotonicity of their difference). The ascending limb of the \(A\)-nullcline and the \(I\)-nullcline are both increasing, giving at most one intersection; the descending limb of the \(A\)-nullcline intersects the still-increasing \(I\)-nullcline at most once (again by monotonicity of the difference). A third intersection can occur in the transition region near the peak. Thus at most three intersections total.
The lower bound uses the intermediate value theorem on the difference \(h(A) = I_\text{null}^A (A) - I^*(A)\). Near \(A = 0\): \(P(0) \approx 0\) (the production machinery requires some ATP), so \(I_\text{null}^A (0) \approx -\beta_0\\/e_a < 0\) while \(I^*(0) = 0\), giving \(h(0) < 0\). As \(A\) increases, \(P(A)\) rises (ascending limb of the unimodal production function) and eventually \(P(A) - \beta_0 > e_a I^*(A)\), making \(h > 0\): the \(A\)-nullcline rises above the \(I\)-nullcline. By the IVT, at least one crossing occurs on the ascending limb (the disease fixed point). On the descending limb at high \(A\), \(P(A)\) falls back below \(\beta_0\) while \(I^*(A) > 0\), so \(h < 0\) again: at least one crossing on the descending limb (the healthy fixed point). This gives at least two crossings. The third crossing (the saddle) arises when bistability inequality holds: the condition ensures that near the peak of the \(A\)-nullcline, the \(A\)-nullcline slope exceeds the \(I\)-nullcline slope, creating a fold where the \(I\)-nullcline overtakes the \(A\)-nullcline from below, producing an additional sign change in \(h\). Numerical evaluation at the stated ME/CFS parameter values (Mathematical Model Details) confirms that three crossings occur. Combined with the upper bound of three: exactly three.
Step 3: Jacobian and eigenvalue classification. At each fixed point \((A^*, I^*)\), the Jacobian is:
$ J = mat(P’(A^*), -e_a; display(), -(k_ + d_a)) mat(f_A, f_I; g_A, g_I) $ {#eq-jacobian-2d}
The trace is \(\tau = f_A + g_I = P'(A^*) - (k_\text{exh} + d_a)\) and the determinant is \(\Delta = f_A g_I - f_I g_A\). Since \(f_I = -e_a < 0\) and \(g_A > 0\) for all \(A^* > 0\), the cross-coupling term \(-f_I g_A = e_a g_A > 0\) always contributes positively to \(\Delta\).
Healthy fixed point \((A_H, I_H)\), on the descending limb of the \(A\)-nullcline: \(P'(A_H) < 0\), so \(f_A < 0\). The trace \(\tau = P'(A_H) - (k_\text{exh} + d_a) < 0\) (both terms negative). The determinant \(\Delta = |P'(A_H)|(k_\text{exh} + d_a) + e_a g_A > 0\). Negative trace and positive determinant: stable (node or spiral).
Disease fixed point \((A_D, I_D)\), on the ascending limb of the \(A\)-nullcline: \(P'(A_D) > 0\), so \(f_A > 0\). The determinant sign follows from the nullcline crossing geometry. At this intersection, the \(I\)-nullcline is steeper than the \(A\)-nullcline (the \(I\)-nullcline crosses from below as \(A\) increases). In terms of Jacobian entries, this means \(g_A \\/ |g_I| > P'(A_D) \\/ e_a\), i.e., \(e_a g_A > P'(A_D)(k_\text{exh} + d_a)\), ensuring \(\Delta > 0\). For the trace: \(\tau = P'(A_D) - (k_\text{exh} + d_a) < 0\) requires \(P'(A_D) < k_\text{exh} + d_a\). Both quantities have units of (time)\(""^{-1}\) in the non-dimensionalised system. This condition holds when \(A_D\) lies in the region where the ascending slope of \(P\) is moderate, which is the case for the stated ME/CFS parameter values (verified numerically in Mathematical Model Details). Negative trace and positive determinant: stable (node or spiral).
Saddle point \((A_s, I_s)\), near the peak of the \(A\)-nullcline: here the \(I\)-nullcline is less steep than the \(A\)-nullcline (the \(A\)-nullcline crosses the \(I\)-nullcline from below as \(A\) increases—the opposite crossing direction). This gives \(P'(A_s) \\/ e_a > g_A \\/ (k_\text{exh} + d_a)\), so \(P'(A_s)(k_\text{exh} + d_a) > e_a g_A\), making \(\Delta = -P'(A_s)(k_\text{exh} + d_a) + e_a g_A < 0\). Negative determinant: eigenvalues have opposite signs, confirming an unstable saddle.
Step 4: Index theory consistency. By the planar fixed-point index theorem applied to the trapping region \(\mathcal{R}\) (Step 1), the sum of indices of all interior fixed points must equal \(+1\). Stable nodes/spirals have index \(+1\); saddles have index \(-1\). Three fixed points with indices \(+1, -1, +1\) sum to \(+1\). \(checkmark\)
The three fixed points are therefore:
- Disease fixed point \((A_D, I_D)\): low ATP, high immune activation. On the ascending limb of the \(A\)-nullcline. The disease attractor.
- Saddle point \((A_s, I_s)\): intermediate ATP, intermediate immune activation. Its stable manifold forms the separatrix—the boundary between the two basins of attraction.
- Healthy-adjacent fixed point \((A_H, I_H)\): high ATP, low immune activation. On the descending limb of the \(A\)-nullcline. A stable attractor, but with a smaller basin of attraction than in the healthy regime.
The disease attractor exists because the positive feedback loop is self-sustaining at ME/CFS parameter values: high immune activation depletes ATP, depleted ATP impairs immune regulation, impaired regulation sustains high immune activation. This cycle does not close in the healthy regime because \(k_\text{recov}\) is sufficiently high to bring activated cells back to resting before ATP is critically depleted.
1.4 Analytical Summary
The reduced 2D system establishes three results for the energy–immune subsystem:
- Both regimes have a healthy attractor: even ME/CFS parameters support a healthy fixed point—but with a smaller basin of attraction. Recovery is theoretically possible without parameter restoration, but requires a large directed perturbation.
- The disease attractor is absent in the healthy regime: at healthy parameter values, only one fixed point exists (The Hill Exponent Is a Modelling Choice). As \(\alpha_\text{CI}\) decreases and \(k_\text{exh}\) increases toward ME/CFS values, the disease and saddle fixed points appear (consistent with a saddle-node bifurcation, though tracing the exact bifurcation curve requires numerical continuation).
- The bistability depends on two structural ingredients: (a) cooperative (\(n \geq 2\)) ATP-dependent immune activation (The Hill Exponent Is a Modelling Choice), and (b) a unimodal production function \(P(A)\) (a consequence of adenylate kinase feedback and PFK-1 product inhibition). Models sharing both ingredients will exhibit qualitatively similar bifurcation structure; models lacking either (e.g., \(n = 1\) Michaelis–Menten activation) will not.
The 2D analysis treats 62 state variables as quasi-steady-state (QSS) parameters. This is justified when the omitted variables equilibrate faster than the two retained variables. For many omitted variables this holds: cytokine half-lives are minutes to hours, membrane potential equilibrates in milliseconds, and metabolite pools turn over in seconds (Eissing et al. 2011). However, some omitted variables—notably epigenetic modifications (\(\tau_\text{epi} \sim\) months) and immune cell exhaustion (\(\tau_\text{exh} \sim\) weeks)—operate on slower timescales than ATP and activated immune cells. These slow variables act as slowly drifting parameters in the reduced system rather than fast slaves; their effect is to slowly move the bifurcation point, not to alter the existence of bistability at a given parameter snapshot. Nonetheless, the qualitative conclusions—bistability exists, healthy parameters support one attractor, ME/CFS parameters support two—should be regarded as analytically demonstrated for the 2D subsystem and conjectured for the full model until verified by numerical continuation analysis (AUTO, MATCONT). Quantitative claims (exact bifurcation thresholds, basin sizes) require the full system.
1.5 Competing Dynamical Hypotheses
The bistable fixed-point model is not the only dynamical explanation for ME/CFS chronicity. Three alternatives deserve explicit consideration.
Continuous parametric shift (no bistability). Under this hypothesis, ME/CFS represents a different region of a single attractor’s basin: parameters (e.g., \(\alpha_\text{CI}\), \(k_\text{exh}\)) drift due to ongoing damage, and symptoms track the parameter shift continuously. Recovery occurs gradually when parameters return. This model predicts that partial treatment should produce partial, graded improvement—and that spontaneous recovery should correlate smoothly with parameter restoration. It is less parsimonious than the bistable model for explaining two clinical observations: (1) many infections of comparable severity resolve without triggering ME/CFS, suggesting a threshold rather than a graded response; and (2) symptoms persist after the apparent trigger resolves, which a purely parametric model explains only by invoking ongoing hidden damage. However, the continuous-shift hypothesis is not decisively falsified by current evidence—a parametric model with rapid parameter change could mimic threshold-like onset, and ongoing subclinical damage (e.g., viral persistence) could explain persistence. Distinguishing the two models requires longitudinal multi-omic trajectories with sufficient temporal resolution to detect whether the state-space path shows a continuous drift or a discontinuous jump.
Limit cycle (oscillatory attractor). Under this hypothesis, the disease state is not a fixed point but a stable limit cycle: the system oscillates between better and worse states, potentially explaining the boom-bust pattern (PEM followed by partial recovery) characteristic of ME/CFS. This is a genuine alternative that the 2D analysis does not exclude. For limit cycles to arise in a 2D system, the disease fixed point would need to be an unstable spiral (positive real part of complex eigenvalues) surrounded by a stable limit cycle (via Hopf bifurcation). Bistability of the Energy–Immune Subsystem shows the disease fixed point is a stable node at the stated parameters; however, for other parameter values a Hopf bifurcation could replace the stable disease node with a limit cycle. The two hypotheses are not mutually exclusive: some patients may sit in a stable disease attractor (steady-state illness) while others orbit a limit cycle (relapsing-remitting course). Distinguishing them clinically requires longitudinal time-series analysis of biomarkers with sufficient temporal resolution to detect oscillatory versus steady-state dynamics.
Epigenetic lock-in (no dynamical bistability needed). T cell exhaustion via epigenetic imprinting is well-documented (Pauken et al. 2016): once T cells commit to an exhausted phenotype through DNA methylation changes, they persist regardless of the original trigger. Under this simpler hypothesis, chronicity reflects epigenetic memory rather than dynamical bistability—no feedback loop is needed, just a one-way epigenetic switch. This hypothesis is compatible with our model (Extended Subsystem Couplings already includes epigenetic consolidation), but would predict that disease persistence depends entirely on the epigenetic state, not on ongoing energy–immune feedback. It is testable: if epigenetic reprogramming (e.g., via demethylating agents) restores immune function without addressing energy metabolism, the epigenetic-only model is supported. If recovery requires simultaneous energy and immune normalisation, the feedback model is supported.
2 Steady-State Multiplicity
The energy–immune coupling alone (energy immune coupling and immune energy feedback) produces bistability for a range of parameter values, as hypothesized in Energy–Immune Coupling. The extended model, with additional positive feedback loops (mast cell–energy, coagulation–oxygenation, epigenetic–parameter), is expected to produce richer attractor structure. Numerical continuation methods (e.g., AUTO, MATCONT) applied to the steady-state equations \(\mathbf{f}(\mathbf{x}^*, \mathbf{\theta}) = \mathbf{0}\) can map the complete bifurcation diagram.
The central prediction of bifurcation analysis, impossible without the mathematical model: the integrated system supports not one but multiple distinct disease attractors, each corresponding to a different ME/CFS phenotype:
- Immune-dominant attractor: high cytokine levels, elevated exhausted T cells, near-normal mitochondrial parameters. Patients in this attractor present with immune-mediated symptoms (lymphadenopathy, sore throat, flu-like malaise) and respond preferentially to immunomodulatory therapy.
- Metabolic-dominant attractor: low ATP, elevated lactate, impaired metabolic flexibility, near-normal immune markers. Patients present with exercise intolerance and PEM as primary complaints and respond to mitochondrial support.
- Neurovascular-dominant attractor: impaired CBF autoregulation, BH₄ depletion, autonomic dysfunction, moderate energy and immune impairment. Patients present with cognitive dysfunction, orthostatic intolerance, and pain as primary symptoms.
- Severe/locked attractor: all subsystems degraded, epigenetic consolidation complete. This attractor has the deepest basin of attraction and the highest intervention threshold for escape.
Each attractor has a characteristic biomarker signature derivable from the model’s steady-state equations. This subtyping is mechanistic: it derives from the mathematical structure of the coupled ODEs, not from statistical clustering of symptoms. The model predicts that patients within the same attractor basin should respond similarly to targeted interventions, providing a rational basis for treatment stratification.
An important caveat: the discrete-attractor subtype model competes with the alternative that ME/CFS heterogeneity forms a continuous spectrum, as proposed by symptom-based severity models (Nacul et al. 2020) (Jason et al. 2005). These are empirically distinguishable: discrete attractors predict clusters in multi-omics data with gaps between them, while a continuous spectrum predicts a smooth distribution. Current data are insufficient to discriminate (most studies lack the multi-domain simultaneous measurements the model requires), and the truth may be intermediate: discrete attractors with noisy parameter variation producing apparent continuity within each basin.
3 Separatrix Topology and Recovery Paths
The separatrices—boundaries between basins of attraction—determine the minimum perturbation required for state transitions. The separatrix between the healthy and each disease attractor defines the “tipping point” for disease onset; the separatrix between a disease attractor and the healthy state defines the recovery threshold. A key model prediction: these separatrices are not symmetric. Due to epigenetic hysteresis (Extended Subsystem Couplings), the recovery separatrix lies further from the disease attractor than the onset separatrix lies from the healthy state. This means recovery requires a larger sustained perturbation (i.e., more intensive intervention) than the original trigger that caused disease onset—formalizing the clinical observation that ME/CFS is easy to trigger and hard to reverse.
The separatrix topology also predicts transitions between disease subtypes: a patient in the immune-dominant attractor who experiences additional mitochondrial insult (e.g., from sustained ROS damage) may transition to the severe/locked attractor without passing through the healthy state. The model maps the conditions under which such inter-attractor transitions occur, identifying the parameter trajectories that lead to progressive worsening.
4 Parameter Sensitivity of Bifurcation Structure
Not all parameters equally influence the bifurcation structure. The codimension of each bifurcation point indicates how many parameters must be simultaneously varied to produce the bifurcation. Low-codimension bifurcations (saddle-node, transcritical) are robust and clinically relevant; high-codimension bifurcations require improbable parameter coincidences. Preliminary analysis of the reduced energy–immune subsystem identifies \(\alpha_\text{CI}\) (Complex I activity) and \(k_\text{exh}\) (immune exhaustion rate) as the two parameters with the strongest influence on the saddle-node bifurcation that creates/destroys the disease attractor. This suggests that interventions targeting these parameters—mitochondrial protectants and immune checkpoint modulators—have the greatest potential to qualitatively change the disease dynamics rather than merely shifting the system within its current attractor.
ME/CFS is not a single disease with a spectrum of severity but a collection of distinct dynamical states (attractors) in a coupled multi-system model. Different attractors correspond to clinically recognizable subtypes with distinct dominant pathophysiology. Disease onset, progression, and recovery are transitions between attractors governed by separatrix geometry. This hypothesis predicts that: (1) patient clustering by multi-omics data should reveal clusters corresponding to model-derived attractors; (2) treatment response should be subtype-specific and predictable from model-derived biomarker signatures; (3) spontaneous recovery occurs preferentially in patients near separatrix boundaries (identifiable by critical slowing down signals, Temporal Evolution and Disease Trajectories); and (4) progressive worsening reflects inter-attractor transitions toward deeper basins. This hypothesis is uniquely generated by mathematical modeling: verbal descriptions of ME/CFS heterogeneity cannot distinguish between a continuous spectrum and discrete attractors, but the mathematical structure of the coupled ODE system makes a specific prediction resolvable by bifurcation analysis.
Certainty: 0.40. Bistability of the energy–immune subsystem is well-supported by the feedback structure; multi-stability of the full 64-variable system requires numerical verification. The correspondence between mathematical attractors and clinical subtypes remains to be validated against multi-omics patient data.