Pathogen diversity emerging from coevolutionary dynamics in interconnected systems

Abstract

The spread of infectious disease and the evolution of antigenically distinct strains are often modeled separately, despite strong feedbacks mediated by host immune memory and heterogeneous contacts. To tackle this challenging problem, we introduce a coevolutionary framework in which transmission occurs on a metapopulation network while mutational exploration of strain space follows a mutation network. In this multiscale model, cross-immunity is encoded by similarity in the latent diffusion geometry of the strain network, so that nearby strains confer partial immune protection. We first identify an effective critical region that controls the transition between extinction, recurrent outbreak episodes, and long-lived endemic persistence, thus characterizing the resulting strain-turnover dynamics. We then derive a replicator-mutator-like equation for strain composition and an explicit dynamical evolutionary landscape induced by the coupling of mutation and transmission. Finally, allowing host heterogeneity to modulate the local mutation structure, we show that spreading across demes can effectively connect otherwise disconnected components of strain space, increasing long-term endemic diversity while producing a non-monotonic change in overall prevalence. Together, our results isolate minimal mechanisms by which immune-mediated competition and network structure can shape antigenic diversification.

Figures (6)

• Figure 1: The evo-SIS framework unfolds over two structures. Locally (Q/K), pathogen strains change according to a mutation network, which also defines cross-infectivity according to the strength of cross-immunity. Globally (EL), the mutation and metapopulation structures become intertwined, as the pathogen spreads across heterogeneous demes and an effective evolutionary landscape emerges.
• Figure 2: (Top) Asymptotic epidemic size $\textbf{(IT}\infty$) and evenness $\textbf{(ET}\infty$) obtained varying $\beta$ in a range comprising the estimated $\beta_\text{crit}^*$, for $\rho \in \{1,2\}$; envelopes represent 10%-90% quantiles. Results are obtained with $M=100$ strains, and are averaged over 200 simulations, in which $\mathbf Q$ (Erdos-Renyi, $p_\text{conn} = 0.12$) is resampled each time. (Bottom) We analyze dynamics for the system \ref{['eq:model_homMF']} in the critical region, $\beta \approx \beta_\text{crit}^* \equiv \left\langle K\right\rangle^{-1}$, where $\langle K \rangle \equiv M^{-1}\sum_{a,b}K_{ab}$, by computing three quantities of interest for different values of cross-immunity strength $\rho$ by sampling 1000 realizations of the same system. (PE) The probability of reaching an endemic steady state at long times quickly saturates to one as $\beta \to \beta_\text{crit}^*$. (DT) Dominance time of a strain is defined as the length of the time interval during which the strain has the largest relative prevalence. The mean dominance time (MDT) is an indicator of the (inverse) rate of turnover between strains. For all values of $\rho$, the MDT decreases sharply as $\beta$ approaches $\beta_\text{crit}^*$ from below. (NO) The mean number of outbreaks peaks as $\beta \to \beta_\text{crit}^*$, and decreases in the super-critical region $\beta > \beta_\text{crit}^*$, as large infectivity quickly drives the system to a steady state. Lines, where present, are mere interpolations, intended to aid visualization. Additional details are reported in Appendix \ref{['app:computational_details']}.
• Figure 3: We compare evo-SIS model \ref{['eq:model_homMF']} predictions (UM, DM, DS) with empirical data on COVID-19 hodcroft2021covariants(ED), where $M=52$ strains are observed; six countries are shown as examples. In all plots, different colors denote the relative prevalence of different strains. We identify likely model parameters by matching mean dominance time and peak epidemic size predicted by the model to the empirical values, calculated globally, which measure respectively $5.42\pm1.69$ weeks and $\left(3.99\pm3.07\right)\times 10^{-2}$ (relative to total population). We report dynamics with a configuration matching these empirical values, both with a directed (DM) and undirected (UM) mutation structure. We also compare the dominant strain trajectory (DS): the undirected case shows resurgence of earlier strains which are not observed in data, while the directed cases --after clearing the earlier strain-- shows a more realistic progression. Additional details are reported in Appendix \ref{['app:computational_details']}.
• Figure 4: Within a single well-mixed deme, in a regime of large transmission rate and weak cross-immunity ($\beta = 21,\rho = 1.4$), the complete dynamics and its approximants give near identical descriptions of relative abundances, and thus of evenness (E). Leading order dynamics fail to predict variations in total incidence (TI), which is however recovered by higher-order approximants. Instead, for a population of heterogeneous demes, we compare eigenvalue density (ED) of the steady state evolutionary landscape \ref{['eq:multiplex_dynamics']}, $\left.\mathbf \Gamma\right|_{t \to +\infty}$, to that of the ensamble of local mutation matrices $\mathbf Q^x$: a spectral gap emerge in the former which is lacking in the latter, indicating the emergence of an interconnected structure. (TO) Strains' evenness increases monotonically with host heterogeneity $N_Q/N$, while epidemic size reaches maximal values for an non-trivial, intermediate value. The plot reports values normalized by their maximum attained value over the full range of $N_Q/N$, in order to display variations on comparable scales. The metapopulation structure is modeled with an Erdos-Renyi network with $p_\text{conn}^A=0.15$, $N =50$, and one connected component. Additional details are reported in Appendix \ref{['app:computational_details']}.
• Figure 5: Deterministic cycles of within-host replicator-mutator dynamics, eq. \ref{['eq:RM_WH_dynamics']}, with an undirected mutation matrix $\mathbf m$. We consider 10 strains, a mutation rate $\varepsilon = 3\times 10^{-3}$, and off-diagonal entries of $\mathbf m$ as equal to $\varepsilon$, multiplied by a random number drawn by an exponential distribution of mean equal to 3. (A) Strain frequencies follow a recurrent pattern, with dominant strains reappearing in cycles (B). At this deterministic level, these are identified (C) by the deterministic mapping in eq. \ref{['eq:within_host_deterministic_successor']}. Dashed lines identify predicted dominant strains both in (B) and (C).