Extinction pathways in dengue-Zika multi-pathogen systems with antibody-dependent enhancement

1 Introduction

Antibody-dependent enhancement (ADE) is an immunological mechanism in which cross-reactive antibodies generated during a primary infection increase susceptibility to, and severity of, subsequent heterologous infections. ADE plays a central role in the epidemiology of dengue virus, where four antigenically distinct serotypes (DENV1–4) co-circulate, and in the interactions between dengue and Zika viruses [1]. A defining feature of ADE is its inherent asymmetry: infection with serotype A may substantially enhance transmission of serotype B, while the reciprocal effect may be weak or absent [2]. This non-reciprocal structure has profound clinical consequences, contributing to severe disease manifestations responsible for approximately 20,000 deaths annually [3].

A fundamental question in infectious disease ecology is: how do multi-pathogen systems with asymmetric interactions transition to extinction through stochastic fluctuations? Understanding this question has direct public health relevance for dengue and Zika control, particularly as climate change expands the geographic range of competent vectors [4]. While deterministic models predict stable pathogen coexistence, real populations occasionally experience local elimination. The mechanisms governing such rare extinction events remain poorly understood, especially when pathogens interact through asymmetric immune-mediated coupling.

Existing theoretical analyses of extinction in multi-pathogen systems face two major limitations. First, they are largely confined to one- or two-pathogen models with symmetric interactions [5], whereas realistic ADE systems are intrinsically high-dimensional and asymmetrically coupled. For a model incorporating four dengue serotypes and Zika virus, the stochastic state space contains five interacting infectious compartments and admits 2⁵−1 = 31 distinct extinction manifolds. Second, dimensional reduction—the standard approach to simplify such systems—can eliminate the very asymmetry and coupling structure that drives extinction dynamics. Consequently, the stochastic mechanisms underlying dengue's remarkable long-term persistence and occasional local elimination have remained unexplored.

In this work, we discover that dengue-Zika extinction occurs through previously unrecognized cascades where elimination of a single “keystone” serotype triggers subsequent pathogen collapses. This cascade phenomenon emerges from the topological structure of the extinction action landscape and reveals why some serotypes are more vulnerable to stochastic elimination than others. We further identify three distinct extinction regimes—sequential, pairwise, and total collapse—with transitions governed by critical values of the basic reproduction number R₀ and mean ADE enhancement η. Counterintuitively, we find that intermediate levels of ADE asymmetry maximize extinction probabilities for specific pathogens, suggesting that immune-mediated coupling can be “tuned” to favor local elimination.

To reach these conclusions, we extend Wentzel-Kramers-Brillouin (WKB) large deviation theory to high-dimensional asymmetrically coupled systems—a necessary technical extension that enables analysis of a previously intractable class of biological problems. This adaptation combines Hamiltonian reduction techniques for asymmetric coupling matrices, analytical characterization of asymmetry-induced bifurcations, and specialized algorithms for computing extinction paths in five or more dimensions. The mathematical framework, while developed for dengue-Zika dynamics, applies broadly to any biological system with asymmetric interactions between components.

Our contributions are threefold. First, we provide the first stochastic explanation for dengue's long-term persistence, showing that increasing serotype diversity lowers the extinction barrier for any single serotype while making global elimination exponentially unlikely (a “safety-in-numbers” effect). Second, we identify keystone serotypes whose targeted elimination could leverage cascade dynamics for more effective disease control. Third, we map conditions under which elimination campaigns are most likely to succeed based on measurable epidemiological parameters. These findings bridge mathematical theory and biological application, offering practical insights for public health planning in dengue-endemic regions.

2 Materials and methods

2.1 Biological background and model formulation

We consider a system of five distinct pathogens: dengue serotypes 1–4 (denoted D₁, D₂, D₃, D₄) and Zika virus (denoted Z). Primary infection follows standard susceptible-infectious-recovered (SIR) dynamics for each pathogen. Secondary infections are subject to antibody-dependent enhancement (ADE), whereby prior immunity to one pathogen increases both susceptibility to and transmissibility of heterologous pathogens. The model includes temporary cross-immunity following heterologous infection, reflecting empirical observations [6]. While dengue and Zika are vector-borne, we simplify to direct transmission for mathematical tractability, a common approach in multi-strain studies [7]. This simplification is a limitation (see Section 8.4) and may affect the quantitative accuracy of extinction barriers, though sensitivity analysis (Appendix D) suggests qualitative features may be robust for some parameter ranges.

2.2 Deterministic SIR-ADE model

The deterministic model tracks the following compartments:

• S: Individuals fully susceptible to all five pathogens.
• I_*i*: Individuals infectious with pathogen i (i = 1, …, 5).
• R_*i*: Individuals recovered from pathogen i and immune to homotypic reinfection.
• M_*ji*: Individuals who have recovered from pathogen j and are currently infected with pathogen i via ADE (j≠i). These individuals are infectious with pathogen i but remain immune to pathogen j; they are distinct from doubly recovered individuals (not explicitly tracked).

The dynamics are governed by Equations 1–4:

dSdt=μN-S∑i=15βiIi/N-μS,(1)

dIidt=βiSIi/N+∑j≠iηjiβiMjiIi/N-(γi+μ)Ii,(2)

dRidt=γiIi-μRi-∑j≠iηijβjRiIj/N,(3)

dMjidt=ηijβjRiIj/N-μMji-ηjiβiMjiIi/N,(4)

where the total population

N=S+∑i=15(Ii+Ri+∑j≠iMij)

is constant. One can verify that summing the right-hand sides of Equations 1–4 yields dN/dt = 0 (all terms cancel pairwise), confirming population conservation. This conservation is a consistency check of the model formulation.

2.2.1 Interpretation of ADE terms

In Equation 2, the term ∑j≠iηjiβiMjiIi/N represents the ADE–mediated force of infection: individuals in class M_*ji* (immune to j, currently infected by i) contribute to the infectious pool for pathogen i at an enhanced rate η_*ji*β_*i*. In Equation 4, the term η_*ij*β_*j*R_*i*I_*j*/N describes the creation of new M_*ji* individuals when a person recovered from i (and thus immune to i) is infected by pathogen j with ADE enhancement η_*ij*. The negative term −η_*ji*β_*i*M_*ji*I_*i*/N accounts for ADE–mediated reinfection of the same individual while still in the M_*ji* state; such individuals remain in M_*ji* (they do not become doubly infected) but contribute to onward transmission. The model does not track double infections beyond the M_*ji* compartment, and there is no double–counting because each individual belongs to exactly one compartment.

2.2.2 Asymmetry measure

The ADE matrix η = [η_*ij*] is generally non–symmetric (η_*ij*≠η_*ji*) and sparse, with enhancement factors typically ranging from 1.5 to 10 [8, 9]. To quantify the degree of asymmetry, we define the scalar measure

α=∥η-ηT∥F∥η∥F,

where ||·||_*F* denotes the Frobenius norm. α = 0 corresponds to a symmetric matrix (η = η^T), and α = 1 corresponds to a perfectly anti–symmetric matrix (η = −η^T). This definition will be used in Section 6.3 to explore the effect of asymmetry on extinction probabilities.

2.3 Stochastic formulation

To study extinction events driven by demographic noise, we construct the continuous-time Markov chain associated with the deterministic system. Let

n=(nS,nI1,…,nI5,nR1,…,nR5,nM11,…,nM55)

denote the state vector, with total population N=∑XnX. The master equation for the probability P(n, t) is Equation 5:

dPdt=∑n′≠n[W(n′→n)P(n′)-W(n→n′)P(n)](5)

where W(n→n′) are the transition rates corresponding to infection, recovery, and demographic events.

For large populations it is convenient to work with the intensive variables x = n/N. Define the rescaled transition rates

w(x′→x)=W(Nx′→Nx)N.

In terms of these, the master equation becomes

∂P(x,t)∂t=N∫dx′[w(x′→x)P(x′,t)-w(x→x′)P(x,t)].(6)

For N≫1 we anticipate that the probability of a large deviation (such as extinction) is exponentially small in N [10]. We therefore insert the WKB ansatz

P(x,t)~exp[-N𝒮(x,t)], x=n/N,(7)

where 𝒮(x, t) is the action function. Substituting Equation 7 into Equation 6 and expanding to leading order in N⁻¹ gives the time-dependent Hamilton–Jacobi equation

∂𝒮∂t+H(x,∇x𝒮)=0,(8)

with the Hamiltonian H(x, p) derived in Appendix A (see Equation 10). For a systematic field–theoretic derivation of the Hamiltonian from the master equation, see Kamenev [11].

Extinction events are exponentially rare in the population size N [12]. On the long timescale over which extinction occurs, the system, conditioned on non-extinction, relaxes to a quasi-stationary distribution. In this regime the action becomes time-independent, ∂𝒮/∂t = 0, and the Hamilton–Jacobi equation reduces to

H(x,∇𝒮)=0.(9)

For a detailed exposition of the WKB method for stochastic population processes we refer the reader to Assaf and Meerson [5, 13]. The explicit form of the Hamiltonian is given by Equation 10 in Appendix A.

2.4 Assumptions

The analysis rests on the following key assumptions:

• Direct transmission (vector dynamics neglected).
• Constant population size (μ balanced).
• Time-independent ADE factors η_*ij*.
• Large population (N≫1) for WKB validity.
• Demographic noise as the only stochastic source.

Their limitations are discussed in Section 8.4.

3 Results

3.1 WKB analysis and Hamiltonian dynamics

3.1.1 Hamilton-Jacobi formulation

Substituting the WKB ansatz (Equation 7) into the master equation and collecting terms of leading order yields the Hamilton-Jacobi Equation 9, where the Hamiltonian is:

H(x,p)=∑iβi×S×Ii(e^{pIi−pS}−1)+ΣiΣ{j≠i}η{ji}βi×M{ ji }×Ii(e^{pIi−pM{ji}}−1)+Σiγi×Ii(e^{pRi−pIi}−1)+μ(1−xS−Σi×Ii−Σi×Ri−ΣiΣ{j≠i}×M{ij})(e^{p s}−1)+Σiμ×Ri(e^{pRi}−1)+ΣiΣ{j≠i}μ×M{ij}(e^{pM{ij}}−1)(10)

where p = ∇𝒮. A detailed derivation is provided in Appendix A.

3.1.2 Extinction boundaries and target manifolds

Extinction of pathogen i corresponds to the hyperplane x_{I_*i*} = 0. More generally, extinction of a subset A⊆{1, …, 5} of pathogens defines an extinction manifold Equation 11:

ℳA={x:xIi=0 ∀i∈A}.(11)

For the five-pathogen system, there are 2⁵−1 = 31 non-empty extinction manifolds. Our objective is to compute the most probable path from the endemic equilibrium x^*, satisfying H(x^*, 0) = 0, to a target extinction manifold ℳ_*A*. The quasi-stationary distribution near x^* implies 𝒮(x^*) = 0, while on ℳ_*A*, 𝒮(x) represents the extinction barrier.

3.1.3 Most probable extinction paths

The most probable extinction paths correspond to zero-energy trajectories of Hamilton's equations:

dxdt=∂H∂p,(12)

dpdt=-∂H∂x,(13)

with boundary conditions:

x(-∞)=x*, p(-∞)=0,(14)

x(+∞)∈ℳA, p(+∞)=0.(15)

These conditions ensure the path starts at the endemic equilibrium (zero momentum, corresponding to deterministic steady state) and ends at the extinction manifold (zero momentum, indicating absorption). The action along such a path is Equations 13–16:

𝒮=∫-∞∞p·dx=∫-∞∞p·dxdtdt,(16)

so that the extinction probability scales as P_ext~exp[−N𝒮].

Solving this high-dimensional boundary value problem requires both analytical reduction (Section 4) and numerical computation (Section 5).

4 Analytical results for special cases

4.1 Two-pathogen subsystem with asymmetric ADE

To isolate the role of asymmetric coupling, we first analyse a reduced system with two interacting pathogens and unidirectional ADE. This setting allows explicit analytical characterization of how asymmetry reshapes extinction pathways while retaining the core biological feature of non-reciprocal immune enhancement.

Proposition 1 (Weak coupling regime). For η−1≪1, the extinction probability for the enhanced pathogen scales as Equation 17:

Pext(2)~exp[-N(𝒮0(2)+(η-1)𝒮1(2)+𝒪((η-1)2))],(17)

where 𝒮0(2) is the extinction action in the uncoupled system (η = 1), and 𝒮1(2)>0. A detailed perturbative analysis confirming 𝒮1(2)>0 is provided in Appendix B.

Remark 1 (Biological interpretation). Proposition 1 demonstrates that weak ADE can increase the extinction barrier for the enhanced pathogen. This counterintuitive effect arises because extinction requires depletion of both directly infected individuals and the ADE-mediated reservoir (M₁₂ compartment). Thus, immune-mediated coupling can paradoxically stabilize pathogens against stochastic extinction.

Proposition 2 (Strong coupling regime). There exists a critical enhancement factor η_*c*>1 such that for η>η_*c*, the most probable extinction path for pathogen 2 involves first extinguishing pathogen 1. The critical value satisfies Equation 18:

𝒮direct(ηc)=𝒮indirect(ηc),(18)

where 𝒮_direct is the action for direct extinction of pathogen 2, and 𝒮_indirect is the cumulative action for sequential extinction (first pathogen 1, then pathogen 2). This bifurcation corresponds to a Maxwell-type transition in the action landscape.

4.1.1 Numerical evaluation of (η_*c*)

For the baseline two–pathogen subsystem (using parameters for DENV4 and DENV2 from Tables 1, 2: β = 0.4day⁻¹, γ = 1/7day⁻¹), we computed 𝒮_direct(η) and 𝒮_indirect(η) numerically. The crossing condition 𝒮_direct(η_*c*) = 𝒮_indirect(η_*c*) yields η_*c*≈2.8. This value marks the transition from the sequential to the pairwise extinction regime and agrees with the phase boundary shown in Figure 1 (vertical dashed line near η≈2.5–3.0).

Table 1. Model parameters and typical ranges.

For Zika virus, transmission and recovery rates are taken from Liu et al. [4]. The baseline ADE matrix is given in Table 1 of the main text.

Table 2. Parameter ranges for dengue serotypes used in numerical simulations.

Remark 2 (Biological interpretation). Proposition 2 shows that sufficiently strong asymmetry forces extinction to proceed via intermediate pathogen loss. Biologically, elimination of one pathogen becomes a prerequisite for eliminating the other, giving rise to cascade extinction. The numerical value η_*c*≈2.8 lies within the empirically observed range for dengue serotypes [2,9], indicating that cascade dynamics are relevant for natural populations.

4.2 Symmetric n-pathogen system

We next consider a symmetric n-pathogen system to disentangle dimensionality effects from asymmetry. Although biological ADE interactions are rarely symmetric, this limit provides analytical insight into how increasing pathogen diversity alters extinction barriers.

Theorem 1 (Invariant manifold). For the symmetric n-pathogen system (β_*i* = β, γ_*i* = γ, η_*ij* = η for all i≠j), the set Equation 19

ℳsym:={ xIi=xI∀i,xRi=xR∀i,xMij=xM∀i≠j }(19)

is an invariant manifold of the Hamiltonian dynamics. On ℳ_sym, the Hamiltonian H reduces to H_sym given by

Hsym=βxSxI(epI-pS-1) +(n-1)ηβxMxI(epI-pM-1) +γxI(epR-pI-1) +μ(1-xS-nxI-nxR-n(n-1)xM)(epS-1) +nμxR(epR-1)

+n(n-1)μxM(epM-1).(20)

Proof. The permutation group S_*n* acts on the pathogen indices by permuting I_*i*, R_*i*, and M_*ij* simultaneously. Because all parameters β_*i*, γ_*i*, η_*ij* (i≠j) are identical, the Hamiltonian H(x, p) is invariant under any permutation σ∈S_*n*. The set ℳ_sym is precisely the fixed–point set of this group action: a point x satisfies x_{I_*i*} = x_{I_*j*}, x_{R_*i*} = x_{R_*j*}, x_{M_*ik*} = x_{M_*jk*} for all i, j, k if and only if it is invariant under all permutations.

To prove invariance, we must show that the Hamiltonian vector field is tangent to ℳ_sym. It suffices to verify that for any two distinct indices i, j,

ẋIi-ẋIj=0, ẋRi-ẋRj=0, ẋMik-ẋMjk=0,

whenever x∈ℳ_sym. Using Hamilton's equations ẋ_*a* = ∂H/∂p_*a* and the explicit form of H (Equation 8 in the main text), we compute these differences. For the I-compartments:

ẋIi=βxSxIiepIi-pS+∑k≠iηkiβxMkixIiepIi-pMki-(γ+μ)xIi.

On ℳ_sym, x_{I_*i*} = x_*I*, x_{M_*ki*} = x_*M* for all k≠i, and η_*ki* = η (symmetric coupling). Moreover, because the momenta p_{M_*ki*} are also equal on the fixed–point set (by the same permutation invariance applied to the conjugate momenta), we have p_{M_*ki*} = p_*M* for all k≠i. Hence

ẋIi=βxSxIepI-pS+(n-1)ηβxMxIepI-pM-(γ+μ)xI,

which is independent of i. Therefore ẋ_{I_*i*}−ẋ_{I_*j*} = 0. Similarly,

ẋRi=γxIiepRi-pIi-μxRi-∑k≠iηikβxRixIkepIk-pRi.

On ℳ_sym, x_{I_*i*} = x_*I*, x_{R_*i*} = x_*R*, p_{R_*i*} = p_*R*, p_{I_*i*} = p_*I*, and η_*ik* = η. The sum contains n−1 identical terms, giving

ẋRi=γxIepR-pI-μxR-(n-1)ηβxRxIepI-pR,

again independent of i. For x_{M_*ij*} (i≠j):

ẋMij=ηijβxRixIjepIj-pRi-μxMij-ηjiβxMijxIjepIj-pMij.

On ℳ_sym, x_{R_*i*} = x_*R*, x_{I_*j*} = x_*I*, x_{M_*ij*} = x_*M*, and all momenta are equal (p_{R_*i*} = p_*R*, p_{I_*j*} = p_*I*, p_{M_*ij*} = p_*M*). Since η_*ij* = η for i≠j, the right–hand side becomes

ẋMij=ηβxRxIepI-pR-μxM-ηβxMxIepI-pM,

which is independent of the pair (i, j). Thus differences vanish. Hence the vector field is tangent to ℳ_sym, proving invariance.

Finally, substituting the equalities x_{I_*i*} = x_*I*, x_{R_*i*} = x_*R*, x_{M_*ij*} = x_*M*, p_{I_*i*} = p_*I*, p_{R_*i*} = p_*R*, p_{M_*ij*} = p_*M* into the full Hamiltonian (Equation 8) and using that there are n choices for i in the I_*i* sums, n for R_*i*, and n(n−1) for M_*ij* terms, we obtain H_sym as given in Equation 8. The x_*S* and p_*S* variables are unchanged.

Remark 3. The reduction to a three–dimensional system (in terms of x_*S*, x_*I*, x_*R*, and x_*M*) is exact on ℳ_sym. This explains why symmetric multi–strain models exhibit degenerate extinction pathways: all pathogens are statistically identical.

Theorem 2 (Extinction probability scaling on the symmetric manifold). For the symmetric n-pathogen system, consider the extinction of a single pathogen (say pathogen 1) starting from the endemic equilibrium. On the invariant manifold ℳ_sym, the extinction action 𝒮_sym(η, n) satisfies, for fixed η>1 and as n → ∞,

𝒮sym(η,n)=C(η)n+D(η)n2+O(n-3),

where C(η)>0 and D(η)>0 are functions defined in Appendix C. Consequently,

∂𝒮sym∂n<0 for all sufficiently largen.

Moreover, the action for simultaneous extinction of all n pathogens satisfies

𝒮total(η,n)=n𝒮sym(η,n)+o(n).

Proof. The proof proceeds by a rigorous asymptotic expansion of the Hamilton–Jacobi equation on ℳ_sym [5, 10]. The reduced Hamiltonian H_sym in Equation 20 depends explicitly on n through the coefficients (n−1), n, and n(n−1). We perform an expansion in 1/n [11]. Define the scaled variables

x̄M=nxM, p̄M=pM,

so that the term (n-1)ηβxMxI(epI-pM-1) becomes ηβx̄MxI(epI-pM-1)+O(1/n) and the demographic terms involving n(n−1)μx_*M* become nμx̄M plus corrections. The full asymptotic analysis is presented in Appendix C. There we show that the Hamilton–Jacobi equation admits a solution of the form

𝒮sym=A(η)n+B(η)n2+O(n-3)

with A(η)>0, B(η)>0. The positivity of A(η) follows from the fact that the leading–order Hamiltonian is a convex function whose minimizer yields a strictly positive action. The positivity of B(η) follows from the second variation of the action being positive definite. The monotonicity ∂𝒮_sym/∂n < 0 is then immediate from differentiation.

For the total extinction action, note that on ℳ_sym the optimal path to eliminate all n pathogens consists of n sequential single–pathogen extinctions. Because the manifold is invariant and the pathogens are indistinguishable, the action for each step is 𝒮_sym(η, n) to leading order. The sub–leading corrections contribute at most o(n). Hence 𝒮_total = n𝒮_sym+o(n). A detailed large–deviation principle argument is given in the Supplementary material [10].

Remark 4. The theorem, as stated, holds for paths that remain on the symmetric manifold. The global minimiser of the action over all paths (including those that break symmetry) is not proven to lie on ℳ_sym; this is an open problem for large n and strong η. Therefore, the biological “safety-in–numbers” interpretation applies strictly to the symmetric subspace. In the fully asymmetric system (as in the dengue–Zika model), symmetry is broken, and the extinction-risk hierarchy is non–degenerate.

5 Numerical methods and results

Rather than merely estimating extinction probabilities, we use numerical action minimization to reconstruct the global geometry of extinction pathways in high–dimensional ADE systems. The numerical results serve to map the topology of the action landscape, identifying basins, saddle points, and barriers that organize extinction dynamics across parameter space.

5.1 Algorithmic pipeline overview

We employ a three–algorithm pipeline specifically designed for high–dimensional, asymmetrically coupled systems:

• Adaptive SVD projection (Algorithm 1) is applied as a pre–processing step to identify the effective low–dimensional subspace containing >95% of the gradient variance. This reduces the search space for subsequent algorithms while preserving asymmetry.
• Minimum action method (MAM, Algorithm 2) computes the minimum–action path and the extinction action value 𝒮 for a given target extinction manifold ℳ_*A* (e.g., extinction of D4 alone). The MAM provides the node weights in Figure 2 and the core values in Table 3.
• Nudged elastic band (NEB, Algorithm 3) locates saddle points on the action landscape that separate two extinction basins. The action barriers Δ𝒮 at these saddles give the edge weights (transition rates) in the pathway network of Figure 2.

Table 3. Parameter ranges for Zika virus used in numerical simulations.

Algorithm 1 Adaptive SVD projection. Steps for dimensionality reduction are listed: Sampling gradients at one thousand points near equilibrium, forming a gradient matrix, computing its thin singular value decomposition, determining effective rank so the ratio of squared singular values exceeds zero point nine five, projecting to the reduced subspace, and solving the reduced problem followed by back-projection. Mathematical notation is included for equations, matrix operations, and transformations.

Algorithm 2 Minimum action method with momentum acceleration. Text-heavy document explaining a computational method for handling local minima and iterative updates, including bullet points, equations, parameter definitions, and step-by-step procedures for initialization, gradient computation, accelerated descent, projection, constraint enforcement, convergence criterion, and output.

Algorithm 3 Nudged elastic band for action saddles. Algorithm instructions presented in five numbered steps for finding saddle points in a dynamical system. Steps outline state initialization, tangent direction definition with two cases, iterative force computations and reparameterization, identification of the saddle point, and saddle verification through Hessian eigenvalue check. Includes detailed formulas and substeps with mathematical symbols for force calculations.

The following subsections detail each algorithm, including convergence criteria and checks for uniqueness.

5.2 Adaptive SVD projection

To handle the full 16–dimensional system efficiently without dimensional reduction artifacts, we use adaptive singular value decomposition.

For the baseline parameters, we found k = 5 retains 96% of the variance, reducing the effective dimensionality from 16 to 5.

5.3 Minimum action method

We solve the Hamilton–Jacobi boundary value problem using the geometric minimum action method [16]. The action functional

𝒮=∫-∞∞[p·x∙-H(x,p)]dt

is discretised on a uniform time grid t_*k* = t_min+kΔt, k = 0, …, K, with Δt = (t_max−t_min)/K.

5.3.1 Grid convergence

For a representative parameter set (baseline parameters, extinction of DENV4 alone), we computed the action 𝒮 with increasing grid resolution K = 50, 100, 200, 500, 1, 000. The action stabilizes to within 10⁻⁸ for K≥500 (see inset in Figure 3, Appendix). All reported results use K = 500.

5.3.2 Uniqueness

For each target manifold, the 50 different initialisations converged to the same action value to within the numerical tolerance (< 10⁻⁸), indicating that the found path is the global minimum for that target. The NEB method (Section 5.4) independently confirmed that no lower–action saddle separates different basins.

5.4 Nudged elastic band for saddle points

To locate saddle points of the action separating extinction basins, we adapt the climbing image nudged elastic band method [17].

The NEB provides the action barriers Δ𝒮 used as edge weights in Figure 2 and in the hysteresis analysis below.

5.5 Parameter specification and implementation

All computations used biologically plausible parameters derived from epidemiological studies (Tables 1, 2). Unless otherwise stated, we set: β_*i* = 0.4 day⁻¹ (range: 0.1–0.6), γ_*i* = 1/7 day⁻¹ (range: 1/5–1/10), μ = 1/(70 × 365) day⁻¹, and N = 10⁶ [14, 15]. The baseline ADE matrix follows the asymmetric structure from Table 1 in the main text.

Algorithms were implemented in Julia 1.9 using DifferentialEquations.jl and Optim.jl. Parallel computation employed 16 threads on an AMD EPYC 7763. The complete codebase is available at https://github.com/richmondbalinia/ADE-extinction.

5.6 Validation and benchmarks

Validation was performed against three benchmarks of increasing complexity.

(i) Single-pathogen SIS model (consistency check on Hamiltonian derivation). The stochastic SIS model has the known WKB extinction action [13]:

𝒮SIS(R0)=ln R0+1R0-1.(21)

We note that the single-pathogen limit of the present model (one pathogen, no ADE, μ>0) is an SIR model with vital dynamics, not an SIS model; the SIR action has no closed form. Equation 21 therefore serves as a check on the correctness of the Hamiltonian derivation and MAM implementation in a related, simpler system. To perform this check, we implemented the stochastic SIS Hamiltonian separately and ran the MAM against the exact values for R₀∈{1.5, 2, 3, 5, 10}. Results are reported in Supplementary Table S4. Maximum relative error was 0.08%, confirming that the MAM implementation is correct for this class of problem.

(ii) Single-pathogen SIR with vital dynamics (direct benchmark for the n = 1 limit of this model). For the SIR-vital dynamics model (the exact n = 1 limit of Equations 1–4 with all ADE terms removed), there is no closed-form WKB action. We therefore validated numerically: for five values of R₀∈{1.5, 2, 3, 5, 10} with γ = 1/7 day⁻¹ and μ = 1/(70 × 365) day⁻¹, we compared the MAM-computed action against extinction probabilities estimated from 5 × 10⁴ Gillespie (direct SSA) simulations at N = 10⁴. The Gillespie estimate of the action is 𝒮^=-ln(P^ext)/N. Results are reported in Supplementary Table S4. The maximum relative difference between MAM and Gillespie estimates was 1.4%, consistent with the expected Monte Carlo sampling error at N = 10⁴.

(iii) Monte Carlo simulations for the five-pathogen ADE model. As reported previously, MAM action values for the full five-pathogen system were compared against Gillespie simulations at N = 10⁴ and N = 10⁵ across 50 parameter combinations, confirming exponential scaling with Pearson r = 0.98 and maximum relative error 0.3% (Supplementary Table S1).

5.6.1 Uncertainty quantification for actions, barriers, and crossings

The global validation against benchmarks establishes the numerical accuracy of the MAM implementation. However, the finer claims in the paper—the precise ranking of pathogen vulnerabilities (Table 3), the barriers Δ𝒮 in the pathway network (Figure 2), and the asymmetry values at which the most endangered pathogen changes (Figure 4)—require additional uncertainty quantification.

To quantify uncertainty on the individual actions 𝒮_*i* in Table 3, we computed 95% confidence intervals from two sources: (i) numerical error from the MAM convergence tolerance (10⁻⁸), which is negligible, and (ii) parameter-induced uncertainty from the ±20% ADE matrix perturbations reported in Section 5.7. For each action 𝒮_*i*, the parameter-induced uncertainty is approximately ±0.05–0.10 (standard deviation across 100 perturbed matrices). Table 3 in the main text should be interpreted with these uncertainties in mind.

Table 4. Extinction actions with 95% confidence intervals derived from the 0.3% maximum relative MAM error (Section 5.6).

CI half-width = 1.96 × 0.003 × 𝒮. No individual-pathogen intervals overlap, confirming the action ranking is statistically supported.

Table 5. Uncertainty analysis for the three Figure 4 crossing points.

𝒮₁, 𝒮₂: action values of the two pathogens at the crossing (equal by definition). δ_CI: 95% CI half-width on the action difference. m: approximate slope of the action difference curve near the crossing (estimated from the MAM curves). δα: resulting uncertainty on the crossing location (Equation 12). Rankings more than δα away from a crossing are statistically supported; rankings within δα of a crossing are ambiguous.

For the barriers Δ𝒮 in Figure 2, the uncertainty propagates from the individual action uncertainties. The cascade gap Δ𝒮 = 8.2 − 6.6 = 1.6 has an approximate 95% confidence interval of [1.38, 1.82]. This remains positive and large relative to the numerical error, so the qualitative conclusion—cascade is exponentially more probable than direct extinction—is robust.

For the Figure 4 crossings, where the identity of the most endangered pathogen changes, the action differences are small. Using the confidence intervals from the ADE matrix perturbations, the 95% confidence intervals for the three crossing points are:

αD4↔D2=0.32±0.05,αD2↔Z=0.57±0.04,αZ↔D1=0.78±0.06.

These intervals reflect the uncertainty in the ADE matrix (Section 5.7) rather than numerical error. The broad qualitative pattern—non-monotonic asymmetry effects with three crossings—is robust, but the precise crossing points should be interpreted with caution. Empirical estimates of η_*ij* with tighter confidence intervals would be required to determine the most endangered pathogen with certainty in any specific setting.

5.7 Numerical results and interpretation

Figure 5 shows the computed action landscape for a three–serotype subsystem. The complex topology reveals: (1) multiple local minima corresponding to distinct extinction states; (2) ridges representing transition states identified by NEB (Algorithm 3); (3) pronounced asymmetry reflecting the non–uniform ADE coupling, which breaks the degeneracy of symmetric models.

Figure 2 presents the complete extinction pathway network. Key features: (1) Cascade extinction – direct extinction is exponentially less probable than the cascade sequence. (2) Keystone pathogens—D4 has the lowest individual action barrier (𝒮 = 1.2); its removal triggers subsequent collapses because it provides strong ADE enhancement to other serotypes. (3) Pathway diversity—alternative sequences exist but are exponentially less probable (e.g., the reverse order has cumulative action 7.9, making it ~e^1.3*N* times less likely).

Figure 1 maps the parameter space into three distinct extinction regimes. The inset shows the hysteresis near η≈3.0: for η < 2.65 the sequential path has lower action; for η>2.95 the pairwise path dominates; between these values both paths are within Δ𝒮 = 0.1, so the dominant mechanism is ambiguous. The hysteresis width is Δη≈0.3. Because empirically observed dengue η values span 1.5–10.0, the system is near or within the ambiguous region, implying that small seasonal or age–dependent variations in ADE strength could switch the dominant extinction pathway.

Figure 4 reveals the counterintuitive relationship between ADE asymmetry and extinction probabilities. Contrary to the expectation that increased asymmetry uniformly raises extinction barriers, D2 and D4 exhibit minima at intermediate asymmetry (α≈0.6). The crossing points indicate that ADE structure can be “tuned” to target specific serotypes for elimination.

The action values in Table 3 and the crossing points in Figure 4 carry numerical uncertainty inherited from the MAM convergence tolerance of 10⁻⁸ and the 0.3% maximum relative error established in Section 5.6. We propagate this uncertainty explicitly.

5.7.1 Uncertainty on the cascade dominance gap

The cascade dominance result rests on Δ𝒮 = 𝒮_direct−𝒮_cum = 8.20 − 6.60 = 1.60. With 0.3% relative error, the 1σ uncertainties are ϵ₁ = 0.003 × 6.60 = 0.020 (on 𝒮_cum) and ϵ₂ = 0.003 × 8.20 = 0.025 (on 𝒮_direct), giving Equation 22

Δ𝒮=1.600±√(0.0202+0.0252)=1.600±0.032.(22)

The lower bound on the gap is Δ𝒮_min = 1.568. The cascade-to-direct probability ratio at this lower bound is exp(1.568N), which for N = 10⁶ has an exponent of 1.568 × 10⁶. The cascade dominance conclusion is therefore safe against the numerical error by many orders of magnitude.

5.7.2 Uncertainty on individual action values

Table 4 reports the 95% confidence intervals on all individual and combined extinction actions, computed as [𝒮(1 − 1.96 × 0.003), 𝒮(1+1.96 × 0.003)]. No two individual-pathogen CIs overlap, confirming that the ordering 𝒮_*D*4<𝒮_*D*2<𝒮_*Z*<𝒮_*D*1<𝒮_*D*3 is statistically supported at the 95% level across the full range of α values examined.

5.7.3 Uncertainty on crossing locations

At each crossing point (α = 0.32, 0.57, 0.78), the two action values are equal by definition, so the values at the crossing are not the issue. The question is how precisely the location of the crossing in α is determined. Let δ𝒮(α) denote the signed difference between the two action curves in question. Near a crossing at α_*c*, δ𝒮(α)≈m(α−α_*c*), where m is the slope of the difference. The 95% CI on δ𝒮 has half-width δCI=1.96ϵ12+ϵ22, where ϵ_*i* are the CI half-widths on the two action values at α_*c*. The uncertainty on the crossing location is then

δα≈δCI|m|.(23)

Table 5 reports these quantities for each of the three crossings in Figure 4.

At each crossing the confidence intervals on the two action values overlap (both by construction and because the curves cross continuously), so the precise value of α_*c* is uncertain by δα≈0.06–0.13. Regions of Figure 4 more than δα away from a crossing have unambiguous pathogen rankings; regions within δα of a crossing should be interpreted as ambiguous with respect to which pathogen is more endangered.

5.7.4 Numerical verification at small n (including n = 4 and 5)

To check how rapidly the asymptotic 1/n approximation becomes accurate, we solved the symmetric endemic equilibrium and the corresponding WKB action for n = 1, 2, 3, 4, 5, 10, 20 using the parameters of Table 2 (β = 0.4, γ = 1/7, μ = 1/(70 × 365), η = 2.5). Table 6 reports the computed 𝒮_sym(n) and the deviation from the leading C(η)/n term. The coefficient C(η) is estimated from the large-n limit as C(2.5)≈0.285.

Table 6. Numerical check of the 1/n scaling for the symmetric n-pathogen system.

𝒮_sym(n) is the WKB extinction action for a single pathogen; C(η)/n is the leading-order prediction with C(2.5)≈0.285. Parameters: β = 0.4, γ = 1/7 day⁻¹, μ = 1/(70 × 365) day⁻¹, η = 2.5.

The table shows that 𝒮_sym(n) decreases monotonically with n. However, for n = 4 and 5, the deviation from the leading C/n term is substantial (approximately 50%–60%). The asymptotic law provides a qualitatively correct description (monotonic decrease), but the quantitative scaling is not purely 1/n until n is much larger than 4. This confirms that the safety-in-numbers effect is present in the symmetric system even for the small n relevant to dengue, but the quantitative scaling is influenced by higher-order terms that are not negligible at n = 4–5.

5.7.5 Conceptual tension with the asymmetric motivation

Theorem 2 is derived for a system with identical parameters across all pathogens—a symmetric idealization that contradicts the very asymmetry we argue is essential for ADE dynamics. We include this result for two reasons. First, it provides the only analytically tractable scaling result for n>2, offering a reference point from which to understand how dimensionality affects extinction barriers. Second, it serves as an upper-bound benchmark: symmetry tends to maximize degeneracy among pathogens and minimize the spread of extinction barriers. In the fully asymmetric system, the degeneracy is broken and the extinction-risk hierarchy becomes non-degenerate, as shown in the numerical results of Sections 5.6–5.7.

Whether a qualitatively similar diversity-persistence relationship holds in the fully asymmetric model is a separate question that requires numerical exploration of the full five-pathogen system. The present paper does not attempt such an exploration; the safety-in-numbers result is presented as a theoretical benchmark, not as a claim about real dengue. The numerical results of the asymmetric system (Table 3, Figure 2) provide the relevant predictions for the dengue-Zika case.

5.8 Sensitivity to ADE matrix uncertainty

The baseline ADE matrix (Table 1) is derived from cross–neutralization and epidemiological studies with inherent measurement uncertainty [1, 8]. To assess the robustness of our qualitative findings, we perturbed the entries of η within plausible error bounds and examined the stability of (i) the keystone pathogen identity, (ii) the dominant extinction cascade sequence, and (iii) the position of Zika virus in the cascade.

We generated 100 random ADE matrices by multiplying each non–diagonal entry η_*ij* by a factor drawn independently from a log–normal distribution with mean 0 and standard deviation σ on the log scale. Two perturbation levels were considered: σ = 0.2 (approximately ±20% variation, as stated in Section 6.4) and σ = 0.4 (±40%). For each perturbed matrix, we recomputed the individual extinction actions 𝒮_*i* (using Algorithm 1) and the dominant cascade sequence (using Algorithm 3 to compute transition barriers). The keystone pathogen was defined as the pathogen with the smallest 𝒮_*i*.

At the σ = 0.2 level, DENV4 remained the keystone pathogen (lowest 𝒮_*i*) in all 100 realizations (100%). The dominant cascade sequence D₄→D₂→Z→D₁→D₃ was preserved in 94 of 100 realizations. In the remaining 6 realizations, the positions of D₂ and Z were swapped (D₄→Z→D₂→D₁→D₃), but D₄ always remained first and D₁, D₃ remained last. Thus, the qualitative cascade structure is robust to ±20% variation in η_*ij*.

At the σ = 0.4 level, DENV4 remained keystone in 94 of 100 realizations. In the remaining 6 realizations, DENV2 became the keystone (lowest 𝒮_*i*). The dominant cascade sequence changed more frequently: the D₄-first sequence occurred in only 68% of realizations; alternative cascades (e.g., D₂→D₄→Z→D₁→D₃ or D₄→Z→D₂→D₁→D₃) were observed. This indicates that the keystone identity is not absolutely fixed when uncertainties are as large as ±40%.

The ADE entries involving Zika (η_*iZ* and η_*Zi*) carry their own uncertainty, with reported ranges of 1.1–2.5 [9]. To examine how this affects Zika's position in the cascade, we performed a separate sensitivity analysis varying only η_*iZ* and η_*Zi* within these bounds (keeping other η_*ij* at baseline). Zika's extinction action 𝒮_*Z* varied by ±0.2 around its baseline value of 1.8, corresponding to a coefficient of variation of approximately 11%. In 85% of realizations, Zika remained the third pathogen in the cascade (after D₄ and D₂). In 10% of realizations, it moved to second position (ahead of D₂), and in 5% it moved to fourth position (after D₁). Thus, Zika's position is moderately sensitive to its cross–ADE parameters, but the qualitative observation that Zika is not the keystone pathogen is robust.

The keystone identity (DENV4) is stable under ±20% perturbations but can change under larger uncertainties (±40%). The cascade sequence is qualitatively robust for ±20% uncertainty. Zika's position in the cascade is moderately sensitive to its own cross–ADE values, but it never becomes the keystone under the ranges examined. These results should be interpreted as a guide to the precision needed in empirical measurements of η_*ij*: a coefficient of variation of ≲20% is sufficient to support the qualitative conclusions of this study.

5.9 Convergence and validation

Figure 3 demonstrates grid convergence. Additional validation is described in detail in Section 5.6 and Supplementary Table S4. Supplementary Table S4 reports MAM values against 𝒮_SIS = ln R₀+1/R₀−1, with maximum relative error 0.08%. Separately, the MAM was validated against Gillespie simulation for the SIR-vital dynamics model (the exact single-pathogen limit of this model), with maximum relative difference 1.4% at N = 10⁴ (Supplementary Table S4). The previously reported sub-percent agreement against the wrong formula is therefore superseded; the validation supports MAM accuracy at the same level.

5.10 Mean extinction times and effective population sizes

Within the WKB framework, the extinction probability scales as P_ext~exp(−N_eff𝒮), and the mean time to extinction (MTE) is given by

τ=1γ+μexp(Neff𝒮),(24)

where 1/(γ+μ)≈7 days is the mean infectious period (the natural timescale of the linearised dynamics near the endemic equilibrium), and N_eff is the effective population size that governs the rare–event statistics. Throughout the paper we have used a census population of N = 10⁶ for illustrative scaling; for DENV4, the smallest individual action barrier is 𝒮 = 1.20, giving N𝒮 = 1.2 × 10⁶, which would imply τ≫10^{500, 000} years—a quantity without physical meaning. This discrepancy arises because the relevant population size for stochastic extinction is not the total human census but the effective size of the dynamically connected infectious pool.

5.10.1 Determinants of N_eff

Several well–established mechanisms reduce N_eff far below the census population in practice:

• Spatial fragmentation. Dengue transmission is spatially heterogeneous; a serotype may circulate in a neighborhood while being absent from adjacent areas [18]. The effective population for local extinction is the size of the locally connected sub–population, not the regional total.
• Seasonal bottlenecks. During inter–epidemic troughs, the infectious count I_*i*(t) for a given serotype can drop to 𝒪(1)–𝒪(10) individuals [1]. At these bottlenecks, the stochastic pool is the number of active infectious individuals, not the susceptible population.
• Between–epidemic fade–out. Empirical records of local serotype disappearance involve periods of zero or near–zero incidence for months to years; such events are initiated by stochastic fade–out from small infectious pools, not from the full endemic equilibrium.

These mechanisms are well recognized in the stochastic epidemic literature [12, 13]. The WKB action 𝒮 retains its predictive value for relative extinction probabilities even when N_eff is small enough that the asymptotic approximation is no longer formally valid, because the exponential form of the MTE is expected to persist, albeit with corrections to the pre-factor and to 𝒮 itself.

5.10.2 MTE as a function of effective population size

Table 7 reports the MTE computed from Equation 24 for DENV4 (𝒮 = 1.20), together with the cascade barriers for DENV2 conditional on D4 extinction (𝒮 = 0.90) and for Zika conditional on D4 and D2 extinction (𝒮 = 0.80). The effective population size N_eff is varied over a range that spans values leading to physically plausible extinction times (years to decades) and values where the WKB asymptotic expansion is formally justified (Neff≥104).

Table 7. Mean time to extinction (MTE) for DENV4, DENV2 (conditional on D4 extinct), and Zika (conditional on D4 and D2 extinct) as a function of effective population size N_eff.

MTEs are computed from Equation 13 with τ₀ = 1/(γ+μ)≈6.998 days. The column “WKB valid” indicates whether Neff≥104, the threshold at which leading–order WKB corrections are 𝒪(N⁻¹) ≤ 10⁻⁴. Cells marked ≫ denote values exceeding the age of the universe (≈1.4 × 10¹⁰ yr).

From Equation 24, the effective population size required to produce a specified MTE τ^* is

Neff=ln(τ/τ0)𝒮.

For DENV4 (𝒮 = 1.20) and target timescales:

τ=1 yr⇒Neff≈3.3,τ=10 yr⇒Neff≈5.2,τ=50 yr⇒Neff≈6.6,τ=100 yr⇒Neff≈7.1.

Thus, the model reproduces observed years–to–decades local extinction timescales only for N_eff≈3–8, depending on the target duration.

5.10.3 Tension between timescale reproduction and asymptotic validity

The numbers above reveal a genuine tension that we state without attempting to resolve:

• The WKB asymptotic expansion is justified for N≫1; we validated the numerical implementation for N = 10⁴ and 10⁵ (Section 5.6), where corrections are 𝒪(10⁻⁴). However, these corrections are no longer small for N_eff~5–10.
• At the same time, N_eff~5–10 is the regime where Equation 24 produces physically plausible extinction timescales.
• The two regimes do not overlap. Consequently, the model does not provide a quantitatively verified prediction of absolute MTE at the effective population sizes required for ecological relevance.

This limitation is inherent to large–deviation theory applied to epidemic extinction [12, 13]. The standard position in the literature is that the WKB action 𝒮 and the associated probability ratios remain informative even when extrapolated outside the formal asymptotic regime, because the exponential form of the MTE is expected to persist at small N, albeit with corrections to 𝒮 and the pre-factor. We do not claim more than this.

5.10.4 Scope of the model's contribution

In light of the above, we clarify what the model establishes and what requires further work:

• Established. The action values 𝒮_*i* and barriers Δ𝒮 computed by the MAM are validated at N = 10⁴–10⁵ with relative error below 0.3% (Supplementary Table S1). The ranking 𝒮_*D*4<𝒮_*D*2<𝒮_*Z*<𝒮_*D*1<𝒮_*D*3 is robust to ±20% perturbations of the ADE matrix (Section 5.8). The probability ratio of cascade to direct extinction, exp[(8.2 − 6.6)N] = exp(1.6N), is correct for any N at which the WKB approximation holds, and the qualitative conclusion—cascade is exponentially more probable than direct extinction—holds for all N>0 because 1.6N>0.
• Requiring further work. Absolute MTE values depend on N_eff, which is not a model output but an empirical quantity reflecting local spatial, seasonal, and demographic structure. Calibrating N_eff for a specific outbreak setting would require either (i) fitting the model to longitudinal serotype–specific incidence data, (ii) Gillespie simulation at small N (where the linear noise approximation of Appendix D3 provides a tractable alternative), or (iii) matching the predicted MTE to observed inter–elimination intervals in a well–characterized surveillance dataset. Until such calibration is performed, the model's contribution is comparative: it establishes the relative ordering of extinction barriers and the cascade topology, not absolute calendar times to extinction.
• Framing in abstract and conclusion. Accordingly, we do not claim that this model explains dengue's long–term persistence in the sense of reproducing observed timescales; rather, it provides a stochastic framework characterizing the relative likelihood of extinction pathways and the structural role of ADE asymmetry in organizing those pathways. The phrase “stochastic explanation for dengue's long–term persistence” has been removed from the abstract and replaced with the more limited claim stated above.

5.11 Hysteresis analysis

The phase boundary near η≈3.0 exhibits hysteresis (inset, Figure 1). Sweeping η upward from 2.5 to 3.5, the sequential path remains the minimum action until η≈2.95; sweeping downward, the pairwise path remains the minimum until η≈2.65. The hysteresis width Δη = 0.3 arises from the Maxwell–type transition described in Proposition 2. For dengue, where η can vary seasonally or by host age [1], this implies that the dominant extinction mechanism may switch dynamically, requiring adaptive control strategies.

Collectively, these numerical findings demonstrate that high–dimensional ADE systems exhibit rich, non–intuitive extinction dynamics. The interplay between dimensionality, asymmetry, and nonlinear coupling creates complex phase spaces with multiple extinction pathways, non–monotonic parameter dependencies, and transitions between qualitatively different dynamical regimes.

6 Biological interpretation of mathematical results

The mathematical analysis in Section 3 yields a set of model–based predictions that, subject to empirical confirmation, may offer biological insights for dengue–Zika systems with ADE. We caution that these interpretations rest on the assumptions of the model (direct transmission, time–averaged ADE, constant population, and large–N asymptotics) and should not be taken as established biological facts without independent field or laboratory validation. The biological implications are summarised in Table 8.

Table 8. Summary of mathematical results and their potential biological implications, assuming the model adequately represents reality.

All implications require empirical testing.

6.1 Keystone pathogens and extinction cascades

The extinction pathway network (Figure 2) suggests that, under our parameter choices, direct extinction is less probable than a cascade through intermediate states. In this specific parameterisation, serotype D₄ has the lowest individual extinction action.

6.1.1 Cautious interpretation

If the empirical ADE matrix resembles our baseline (Table 1), then D₄ might act as a keystone pathogen. Its removal could, in principle, trigger subsequent declines in D₂ and Zika. However, this conclusion is sensitive to the uncertain values of η_4*j* (see Section 5.8) and to the direct–transmission approximation. Field studies measuring cross–enhancement from D₄ to other serotypes are needed before any operational recommendation can be made.

6.2 Extinction regimes and control strategies

The phase diagram (Figure 1) distinguishes three regimes based on the magnitude of ADE enhancement η. Within the model, the dominant extinction pathway changes with η.

6.2.1 Cautious interpretation

If local ADE parameters fall into the sequential regime (η < 2.5), then targeting the most vulnerable pathogen individually may be sufficient. If they lie in the pairwise regime (2.5 < η < 4.0), coordinated control of enhancing pairs might be necessary. For very high η (>4.0), mass vaccination could become relatively more effective. However, empirical estimates of η have substantial uncertainty (typical ranges 1.5–10.0), and the boundaries themselves depend on other parameters (e.g., β, γ). Therefore, these regimes should be viewed as qualitative guides rather than precise thresholds for decision–making.

6.3 Asymmetry as a tunable parameter

Figure 4 shows a non–monotonic relationship between the asymmetry measure α and extinction actions for D₂ and D₄ within the model.

6.3.1 Cautious interpretation

This result suggests that, in principle, modulating ADE asymmetry (rather than eliminating it) might increase elimination probabilities for some serotypes. The model predicts an optimal asymmetry around α≈0.6. However, no current intervention is known to achieve such tuning. Moreover, the model does not account for the increased disease severity associated with ADE; any attempt to modulate asymmetry would need to balance elimination benefits against clinical risks. This should be considered a hypothesis for future research, not an actionable recommendation.

6.4 Safety in numbers: a possible explanation for dengue persistence

Theorem 2 demonstrates that on the symmetric invariant manifold, the single–serotype extinction action decreases as 1/n, while the total extinction action increases with n.

6.4.1 Cautious interpretation

This provides a possible stochastic mechanism for the observation that multi–serotype dengue systems often persist globally despite local serotype extinctions. However, the theorem relies on the symmetric manifold assumption; real dengue systems are highly asymmetric. Whether a similar scaling holds in the full asymmetric model is unknown. Moreover, the result does not account for serotype replacement or evolution [21]. Hence it should be seen as a theoretical insight rather than a definitive explanation.

7 Potential public health implications

Based on the mathematical analysis, we outline several model–derived hypotheses that, if confirmed by empirical studies, might eventually inform public health strategies. All recommendations are explicitly conditional on model assumptions and should not be implemented without independent validation. The model-derived hypotheses are summarised in Table 9.

Table 9. Model–derived hypotheses with potential public health relevance (not established guidelines).

7.1 Illustrative (not prescriptive) three–phase strategy

As an illustration of how the cascade dynamics might inform thinking, we sketch a hypothetical three–phase approach. This is not a recommendation for implementation but a scenario to test in future modeling and pilot studies.

• Phase 1 (Keystone targeting): In regions where serotype D₄ is suspected to be a keystone (based on local ADE measurements), intensified surveillance and targeted vaccination against D₄ could be considered as a research priority.
• Phase 2 (Cascade monitoring): If D₄ suppression is achieved, monitor D₂ and Zika incidence for evidence of secondary declines, as predicted by the model.
• Phase 3 (Mop–up): Should serotype diversity fall to 2–3 serotypes, population–wide vaccination might become more feasible, but only after careful cost–effectiveness analysis.

7.2 Critical caveats and unresolved questions

The following limitations mean that any public health inference must be highly provisional:

• ADE matrix uncertainty: The empirical values of η_*ij* have wide confidence intervals (see Section 5.8). Keystone identity may change with different estimates.
• Direct transmission approximation: The model abstracts vector dynamics. A full two–host model could alter extinction barriers and cascade sequences.
• Time–averaged ADE: Enhancement factors vary with time since infection and host age; this may introduce memory effects not captured here.
• Serotype replacement: The model does not account for evolution of new serotypes or replacement following elimination, which could reset the ADE matrix.
• Severe disease risk: Our analysis focuses only on extinction probability, not on clinical outcomes. Interventions that increase extinction might also increase severe dengue risk if ADE is inadvertently enhanced.

7.3 Need for empirical validation

All model predictions—keystone identity, regime boundaries, optimal asymmetry, and safety–in–numbers scaling—require empirical testing. We recommend:

• Longitudinal serological studies to estimate η_*ij* with tighter confidence intervals, particularly for D₄ as a source of enhancement.
• Controlled intervention trials (e.g., targeted D₄ vaccination in well–characterized populations) to test cascade predictions.
• Spatially explicit studies to examine whether local serotype extinctions follow the predicted sequence.

Until such evidence is available, the mathematical results should be viewed as hypotheses to guide further research, not as operational guidelines for public health practice.

8 Discussion

Much of the existing literature on stochastic epidemic extinction using WKB and large deviation methods has focused on low-dimensional systems, typically involving one or two pathogens and a small number of absorbing states [5, 13]. In these settings, extinction is often governed by a single dominant optimal path, and rare-event dynamics can be understood in terms of perturbations around a unique extinction boundary. These assumptions have been successful but implicitly rely on low dimensionality, weak or symmetric coupling, and limited competition between extinction outcomes.

The ADE system studied here exhibits a network of competing extinction pathways separated by non–trivial action barriers. This multiplicity arises from the combination of high dimensionality and asymmetric coupling. Rather than a single optimal path, extinction dynamics in the model are organized by a landscape of locally optimal trajectories connecting multiple absorbing manifolds. As a result, extinction proceeds through cascades of intermediate states, and small parameter changes can induce abrupt switches in the dominant extinction mechanism within the model. Whether these features occur in real biological systems is an open question; they represent model predictions that require empirical testing.

8.1 Mathematical implications (within the model)

Our work adapts and extends existing WKB large deviation theory to a class of biological problems that has received limited attention: high–dimensional asymmetrically coupled multi–pathogen systems. The following observations are strictly mathematical properties of the model; their biological relevance is not established.

First, we have demonstrated that WKB methods, traditionally applied to low–dimensional systems, can be extended to analyze extinction in multi–pathogen systems with asymmetric coupling, at least for the moderate dimensions (up to five pathogens) considered here. The Hamiltonian formulation provides a framework for computing extinction probabilities for systems with interaction structures that include ADE, and in principle could be adapted to other forms of pathogen interaction (cross–immunity, resource competition, etc.), though each adaptation would require its own validation.

Second, within the model, high dimensionality leads to a more complex action landscape compared to low–dimensional approximations. The landscape develops multiple basins and saddle points corresponding to distinct extinction pathways (Figure 5). This complexity suggests that mean–field reductions or pairwise approximations may not capture all features of extinction dynamics in multi–pathogen systems, though the degree of approximation error has not been systematically quantified. The cascade behavior observed in Figure 2—where extinction proceeds through intermediate states—has no direct analog in one– or two–pathogen versions of the model.

Third, the asymmetric nature of ADE coupling breaks the symmetry assumed in many existing rare–event analyses. While fully symmetric systems exhibit degenerate extinction pathways (Theorem 1), asymmetry lifts this degeneracy, creating a hierarchy of extinction probabilities. The non–monotonic relationship between asymmetry and extinction probability (Figure 4) is a mathematical property of the model; whether it has any counterpart in real biological systems is unknown.

8.2 Potential biological insights (model predictions)

From a biological perspective, the model generates several predictions that could, in principle, be tested empirically. We emphasize that these are predictions of the model under its stated assumptions, not established facts.

• Keystone pathogens: Under the baseline ADE parameterization (Table 1), the model predicts that serotype D₄ would be the first to go extinct stochastically, and its extinction would increase the probability of subsequent extinctions of D₂ and Zika. This prediction depends sensitively on the ADE matrix; if empirical estimates of η_4*j* differ substantially, the keystone identity could change.
• Extinction regimes: The model predicts three qualitatively different regimes (sequential, pairwise, total collapse) separated by critical η values. For the parameter values typical of dengue, the model places the system in or near the pairwise regime, suggesting that coordinated fluctuations may be important. Whether these regime boundaries correspond to observable phenomena in real populations is not known.
• Asymmetry tuning: The model predicts that intermediate levels of ADE asymmetry (α≈0.6) maximize extinction probabilities for D₂ and D₄. This is a purely theoretical observation; no intervention currently exists to modulate ADE asymmetry, and even if it did, the trade–off with disease severity would need careful evaluation.
• Safety in numbers: On the symmetric invariant manifold, the model exhibits a scaling where single–pathogen extinction action decreases with the number of serotypes n, while total extinction action increases. This offers one possible stochastic explanation for the persistence of multi–serotype systems, but the assumption of symmetry is strong; the effect may be different in the full asymmetric model.

8.3 Comparison with existing literature

Our model extends the deterministic coexistence results of Ferguson et al. [19], who showed that ADE promotes multi–strain persistence, by examining stochastic mechanisms through which this persistence might eventually be broken. Bianco et al. [20] demonstrated cross–immunity and ADE interplay in two–strain systems; our work extends this to five pathogens with asymmetric coupling and provides extinction probability estimates rather than qualitative dynamics. Wearing and Rohani [18] identified ecological drivers of dengue epidemics; our framework provides a stochastic complement, quantifying how these drivers might translate to extinction risk under the model's assumptions. None of these prior works addressed cascade extinction or keystone pathogen concepts in stochastic multi–strain systems. However, direct comparison is complicated by differences in model structure (e.g., direct vs. vector transmission, time–averaged vs. time–varying ADE).

8.4 Limitations and future extensions

Several limitations of the current analysis should be acknowledged. These limitations mean that the model predictions should be treated with caution until further validated.

• Direct transmission assumption: Dengue and Zika are vector–borne. Incorporating explicit mosquito dynamics would introduce additional dimensions and timescales, and would change the stochastic noise structure in ways that a single effective rate cannot fully capture. Sensitivity analysis (Appendix D) suggests that this approximation may change action values by less than 15% for some parameters, but this bound does not guarantee that qualitative cascade structures remain unchanged.
• Time–averaged ADE: Enhancement factors η_*ij* vary with time since primary infection and host age [1]. Time–dependent ADE could introduce memory effects that shift keystone identities and alter cascade sequences. Age–structured modeling would be required to assess this.
• ADE matrix uncertainty: The empirical ADE matrix (Table 1) is derived from cross–neutralization and epidemiological studies with inherent measurement uncertainty. Sensitivity analysis (Section 5.8) shows that ±20% variation in η values shifts action values by less than 8% in the model, but qualitative changes (e.g., keystone identity) can occur at larger perturbations (±40%). More precise empirical estimates of η_*ij* are needed.
• Leading–order WKB: The current analysis uses leading–order WKB asymptotics. Next–order corrections [13] could refine probability estimates near critical points (e.g., the hysteresis region near η≈3.0) where the leading approximation may be less accurate.
• Serotype replacement: The model does not account for evolution of new serotype variants following elimination, which could reset the ADE matrix and alter post–elimination dynamics.
• Population size: The WKB approximation assumes N≫1. While we have verified results for N = 10⁴ and N = 10⁵, extrapolation to very small populations (e.g., isolated communities) would require different methods.

8.5 Mathematical approach: adaptation, not invention

The mathematical framework developed here adapts and extends existing WKB large deviation theory to a previously underexplored class of biological problems—high–dimensional asymmetrically coupled multi–pathogen systems. Rather than claiming fundamentally new mathematics, we emphasize three methodological contributions that enabled the model–based predictions above:

• Extension to asymmetric coupling matrices: Classical WKB analyses in epidemiology often assume symmetric or weak coupling. We develop Hamiltonian reduction techniques that preserve asymmetry structure, revealing bifurcations and non–monotonic effects that would be invisible under symmetric approximations.
• Multiple competing extinction manifolds: Standard analyses typically consider a single absorbing state. Our framework accommodates 2^*n*−1 competing extinction targets, enabling identification of cascade dynamics and candidate keystone pathogens within the model.
• Numerical methods for high–dimensional action minimization: We adapt the geometric minimum action method and nudged elastic band algorithms to systems with five or more dimensions and complex manifold structures, providing computational tools for problems that have been computationally challenging.

These extensions, while grounded in existing theory, represent necessary adaptations for the specific features of ADE systems. The framework may serve as a template for analyzing other high–dimensional biological systems with asymmetric interactions, but each such application would require its own validation.

8.6 Need for empirical validation

All model predictions—keystone identity, regime boundaries, optimal asymmetry, safety–in–numbers scaling—require empirical testing. We recommend:

• Longitudinal serological studies to estimate η_*ij* with tighter confidence intervals, particularly for η_4*j* (DENV4 as source of enhancement).
• Controlled intervention trials (e.g., targeted D₄ vaccination in well–characterized populations) to test cascade predictions.
• Spatially explicit studies to examine whether local serotype extinctions follow the predicted sequence.
• Comparative studies across regions with different serotype compositions to test the safety–in–numbers scaling.

Until such evidence is available, the mathematical results should be viewed as hypotheses to guide further research, not as operational guidelines for public health practice.

9 Conclusion

We have developed a theoretical framework for analyzing rare extinction events in a multi–pathogen system with antibody–dependent enhancement, using dengue and Zika as a case study. By combining WKB asymptotics with Hamiltonian dynamics and three numerical algorithms, we have derived model–based predictions for extinction probabilities and pathways in systems of up to five interacting pathogens. The principal model predictions, subject to the limitations discussed above, are:

• Cascade extinction: Within the model, extinction typically proceeds through cascades of intermediate states rather than direct transitions to complete elimination. Under the baseline parameterization, serotype D₄ appears as a candidate keystone pathogen.
• Extinction regimes: The model exhibits three distinct regimes (sequential, pairwise, total collapse) separated by critical values of the ADE enhancement factor η. Most empirically observed dengue parameters fall in or near the pairwise regime.
• Asymmetry effects: The model predicts a non–monotonic relationship between ADE asymmetry and extinction probability, with intermediate asymmetry (α≈0.6) maximizing extinction for certain serotypes.
• Safety-in-numbers scaling (on the symmetric manifold): In a symmetric n-pathogen system, the individual extinction action decreases as 1/n while the action for global elimination grows as n·𝒮_sym. This provides a possible mechanistic reference point for understanding how serotype diversity could simultaneously promote local vulnerability and global persistence. Whether a similar relationship holds in the fully asymmetric dengue-Zika system is an open question; the result should therefore be regarded as a theoretical insight derived from a symmetric idealization, not as a direct explanation of empirical dengue dynamics.

These predictions have potential implications for disease control, but only if they survive empirical validation. The identification of candidate keystone pathogens suggests hypotheses for targeted intervention strategies; the phase diagram offers hypotheses about when different control approaches might be most effective; and the scaling relations generate hypotheses about how population size and serotype diversity interact to determine extinction risk.

Looking forward, the framework can be adapted to analyse other high–dimensional biological systems with asymmetric interactions, such as competing ecological species or within–host pathogen dynamics. However, each such adaptation would require its own validation and careful consideration of model assumptions.

In summary, this work bridges mathematical theory and biological application, providing model–generated hypotheses about rare events in complex systems and computational tools for exploring pathogen coexistence and extinction. As climate change and globalization continue to alter disease landscapes, such theoretical frameworks—when properly validated—may become useful for understanding multi–pathogen threats. The strong dependence of our predictions on uncertain parameters (especially the ADE matrix) underscores the need for continued empirical measurement and model refinement.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2026.1841634/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Feng Rao, Nanjing Tech University, China

Reviewed by: Hsuan-wei Wayne Lee, Institute of Sociology, Academia Sinica, Taiwan

Ifeanyi Onah, University Health Network (UHN), Canada