# How Host Mobility Patterns Shape Antigenic Escape During Viral-Immune Coevolution

Natalie Blot, Caelan Brooks, Daniel Swartz, Eslam Abdelaleem, Martin Garic, Andrea Iglesias-Ramas, Michael Pasek, Thierry Mora, Aleksandra Walczak

► To cite this version:

Natalie Blot, Caelan Brooks, Daniel Swartz, Eslam Abdelaleem, Martin Garic, et al. How Host Mobility Patterns Shape Antigenic Escape During Viral-Immune Coevolution. PRX Life, 2025, 3 (2), pp.023011. 10.1103/PRXLife.3.023011. hal-05353331

HAL Id: hal-05353331

https://cnrs.hal.science/hal-05353331v1

Submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
How Host Mobility Patterns Shape Antigenic Escape During Viral-Immune Coevolution

Natalie Blot,<fnref n="1" />,* Caelan Brooks,<fnref n="2" />,*, Daniel W. Swartz,<fnref n="3" />,*, Eslam Abdelaleem,<fnref n="4" /> Martin Garic,<fnref n="5" />

Andrea Iglesias-Ramas,<fnref n="6" /> Michael Pasek,<fnref n="4" />,<fnref n="7" /> Thierry Mora,<fnref n="8" />,† and Aleksandra M. Walczak,<fnref n="8" />,‡

<figure><img src="image_2.png" /></figure>

(Received 22 October 2024; accepted 12 May 2025; published 3 June 2025)

Viruses like influenza have long coevolved with host immune systems, gradually shaping the evolutionary trajectory of these pathogens. Host immune systems develop immunity against circulating strains, which in turn avoid extinction by exploiting antigenic escape mutations that render new strains immune to existing antibodies in the host population. Infected hosts are also mobile, which can spread the virus to regions without developed host immunity, offering additional reservoirs for viral growth. While the effects of migration on viral persistence have been investigated, we know little about how antigenic escape coupled with migration changes the survival and spread of emerging viruses. By considering the two processes on equal footing, we show that on short timescales an intermediate host mobility rate increases the survival probability of the virus through antigenic escape. We show that more strongly connected migratory networks decrease the survival probability of the virus. Using data from high traffic airports, we argue that current human migration rates are beneficial for viral survival.

DOI: 10.1103/PRXLife.3.023011

I. INTRODUCTION

Viruses and their hosts have coevolved since the earliest
form of cellular life [1]. In humans and other jawed ver-
tebrates, infected hosts produce a specific adaptive immune
response through the activation and proliferation of specific B
and T cells that prevent viral spreading and kill infected cells.
Once the infection has been cleared, hosts retain an immune
memory in the form of memory B and T cells, enabling rapid
response in the case of reinfection [2,3]. This long-lasting
immune protection places strong selective pressures on cir-
culating pathogenic strains, driving the coevolution of viral
surface proteins in the face of ever-adapting immune protec-
tions [4–8].

The evolution of the influenza virus is an example of
the coevolutionary feedback experienced by viruses and the
population of its host immune systems [8–11]. Globally cir-
culating influenza strains compete for hosts, which produces

a selective force acting on the genetic variation within the
influenza population [12,13]. Together, selection and viral
mutations occurring at high rates fuel the rapid emergence
of new antigenic variants which in turn lead to the evolution
of host immune systems that acquire protection [5,8,14,15].
Recent in-laboratory expression of historical 20th century
influenza strains shows that the virus did undergo antigenic
escape [15,16].

The influenza virus has been stably evolving in humans for
a long time [10], finding itself now in a coevolutionary steady
state [17]. While selection pressures acting on influenza B
strains have led to the coexistence of two stable lineages,
Victoria and Yamagata, the evolution of influenza A/H3N2
has resulted in one stable lineage since the late 1960s. Dimen-
sionality reduction of strains within influenza A/H3N2 [18,19]
shows that the antigenic escape is well described in terms of a
one-dimensional traveling antigenic wave [19–24].

Antigenic wave descriptions of coevolution assume a
single well-mixed host population, ignoring spatial effects.
Large-scale spatial models of epidemiology have been instru-
mental in guiding policies and explaining epidemic dynamics
[25–28], but usually do not include antigenic escape. Host
migration is known to be an important factor in influenza
spreading [12,29–32]. The effect of fragmenting a popula-
tion into many subpopulations is known to strongly affect
the persistence of species in different ecological scenarios
[33,34], and a similar fragmentation of potential hosts can

<footnote n="1">¹Gulliver Laboratory, UMR No. 7083, CNRS, ESPCI Paris, Université PSL, 75005 Paris, France</footnote>

<footnote n="2">²Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA</footnote>

<footnote n="3">³Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA</footnote>

<footnote n="4">⁴Department of Physics, Emory University, Atlanta, Georgia 30322, USA</footnote>

<footnote n="5">⁵Laboratoire Jean-Perrin, Institut de Biologie Paris-Seine, CNRS, Sorbonne Université, 75005 Paris, France</footnote>

<footnote n="6">⁶Laboratoire Physique des Cellules et Cancer, UMR No. 168, CNRS, Institut Curie, PSL Research University, Sorbonne Université, 75005 Paris, France</footnote>

<footnote n="7">⁷Initiative in Theory and Modeling of Living Systems, Emory University, Atlanta, Georgia 30322, USA</footnote>

<footnote n="8">⁸Laboratoire de Physique, École Normale Supérieure, CNRS, PSL University, Sorbonne Université, and Université Paris Cité, 75005 Paris, France</footnote>

<footnote>*These authors contributed equally to this work.</footnote>

<footnote>†Contact author: thierry.mora@phys.ens.fr</footnote>

<footnote>‡Contact author: aleksandra.walczak@phys.ens.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematic illustrating the different spaces encompassed by the model. (a) Network space. Demes of constant population sizes communicate through host migration. (b) Antigenic space. Within each deme, the dynamics follow the model described in Eqs. (1) and (2). (c) Geographical space. The connectivity of demes within the network is facilitated by the movement of infected people in geographic space. In Sec. III E we connect the spaces in (a) and (b) to the geographic space described in (c).</figcaption></figure>

change epidemiological features of a disease, such as the epidemic threshold [35]. However, few models capture all three features—epidemiological dynamics, antigenic escape, and the spatial structure of populations—explicitly. References [36,37] treat epidemiological dynamics and antigenic escape together in a well-mixed population, while Ref. [38] explores the evolutionary impact of punctuated antigenic escape. Others treat population structure and epidemic dynamics [35,39] or network structure and evolutionary impact without a host-disease model [40,41]. Ecologists have studied “ecoevo coupling” with great interest in the past few decades as mounting evidence shows the importance of such phenomena, but often in well-mixed populations and at steady state. Finally, spatial models of stochastic evolution have coupled small population sizes with spatial structure, allowing individuals to migrate between spatially distinct patches which we call demes, and explored the role of migration on population survival [42–45], but these models have not been applied to coevolutionary settings.

Here we focus on antigenic escape in a stochastic coevolutionary model to study the effects of host migration on the survival of a new viral strain following its outbreak. We ask how the interplay between host spatial migration and evolution within an abstract antigenic space influences the survival probability of the virus (Fig. 1).

## II. MODEL

We describe viral-immune coevolution in terms of a stochastic model of antigenic drift coupled to epidemiological dynamics [24] in a population of identical hosts structured in demes. We focus on mutations that are neutral besides their interaction with the immune system and do not affect life-history traits like transmission or recovery rates. Nevertheless, these mutations can reduce the binding affinity between

antibodies and antigens, effectively allowing the evolved strains to escape host immunity.

Individuals are well-mixed within each deme, and demes are coupled by migration between them [Fig. 1(a)]. We describe the interactions between viruses and hosts with a susceptible-infected-recovered (SIR) model. To effectively describe the space of all possible antigenic strains, we introduce a one-dimensional antigenic space [46] and label the viral strains and host immunity by a continuous antigenic coordinate x. While antigenic space is generally of higher dimension, canalized evolution makes it effectively one dimensional [18,19] if we ignore rare splitting events [22,23]. The density of hosts infected by strain x at a given time t in deme i is $n_i(x, t)$, and the fraction of hosts in deme i who are susceptible to x is denoted by $S_i(x)$. Upon introduction into a deme that is entirely susceptible to strain x [$S_i(x, t) = 1$], the virus will grow exponentially with a characteristic transmission rate $\beta$. An infected host will mount an adaptive immune response and clear the infection with a recovery rate $\gamma$ or die with rate $\alpha$ such that $\partial n_i(x, t)/\partial t = [\beta S_i(x, t) - \gamma - \alpha]n_i(x, t)$. For influenza infections in otherwise healthy human hosts, recovery is the most common outcome, so we set $\alpha = 0$.

A viral strain x accumulates many antigenic mutations that are small in magnitude such that it effectively diffuses through the antigenic space with a diffusion coefficient $D = (\delta x)^2/2\delta t$, where $\delta x$ is the typical effect size of a single mutation occurring on a typical timescale $\delta t$. The viral population density in antigenic space changes due to fitness, migration, and mutations,

$$ \frac{\partial n_i(x, t)}{\partial t} = [\beta S_i(x, t) - \gamma - \alpha]n_i + \sum_{j \neq i} (K_{ij}n_j - K_{ji}n_i) + D \frac{\partial^2 n_i}{\partial x^2} + \sigma \sqrt{n_i} \xi_i(t), \quad (1) $$
where $K_{ij}$ is the migration rate from deme $j$ to deme $i$ and $\xi_i(t)$ is a Gaussian white noise of unit amplitude. Unless specified otherwise, $K_{ij}$ is symmetric so that all demes have equal size at steady state. The parameter $\sigma$ corresponds to the standard deviation of the reproductive number (see the Appendix E).

In response to the viral strains, hosts evolve a population-level immune density in each deme $i$, $h_i(x, t)$. For numerical simplicity, we assume that each deme has a constant number of hosts, $N_h$. Each host has $M$ immune protections and a host infected by an unrecognized strain $x$ randomly replaces an existing protection with immunity to $x$, leading to the dynamics

$$ \partial_t h_i(x, t) = \frac{1}{MN_h} [n_i(x, t) - N_i(t)h_i(x, t)], \quad (2) $$

where $N_i(t) = \int n_i(x, t)dx$ is the total number of infected hosts in deme $i$. We assume immunity decays exponentially with antigenic distance, with characteristic length $r_0$. The population-level immune coverage is defined as

$$ c_i(x, t) = \int h_i(y, t)e^{-|x-y|/r_0} dy. \quad (3) $$

The probability for a potential host to be susceptible to $x$ is then given by the probability that none of their $M$ memories protect them from strain $x$,

$$ S_i(x, t) = [1 - c_i(x, t)]^M. \quad (4) $$

This model was previously studied at steady state in the case of a single deme [24]. It was shown to converge to a stationary solution in a frame moving with constant speed. The speed, size, and width of the wave depends on the cross-reactivity range of the immune protection. Two tractable regimes emerge when two timescales of the problem are well separated: the typical escape time of a mutating virus and the basic doubling time of the virus in the absence of immunity, $\kappa = r_0^2(\beta - \alpha - \gamma)/D$. For relatively small ratios of these timescales ($\kappa \ll 10^3$), the coevolutionary dynamics is described by a Fisher-Kolmogorov-Petrovsky-Piskunov wave, and for large ratios ($\kappa \gg 10^3$) it is described by a linear-fitness wave. Our simulations were done close to the crossover regime, with parameter choices resulting in $\kappa \sim 10^3$.

In this paper we want to study the outbreak of a new strain to which the population is entirely susceptible. To do so, we initiate all the simulations with no immune coverage $h_i(x, t = 0) = 0$ and with a small number of infected hosts appearing in a single deme, $n_1(x, t = 0) > 0$ and $n_{i \neq 1}(x, t = 0) = 0$.

### III. RESULTS

#### A. One deme

To build intuition, we first study the simplest model that includes a single deme, eliminating spatial structure entirely. To study the effect of antigenic mutations on the outbreak dynamics, we numerically integrate Eqs. (1) and (2). Our initial antigenic distribution is isogenic with $N_0 = 100$ infected individuals, all carrying the ancestor strain $x = 0$, $n(x, t = 0) = N_0\delta(x)$. We find that the survival of the viral population

is highly sensitive to noise stemming from population number fluctuations [Fig. 2(a)]. The total number of infected individuals

$$ N_1(t) = \int n_1(x, t)dx \quad (5) $$

follows two possible scenarios during an initial outbreak. After quickly rising to an outbreak peak, the infected host count either falls to zero as hosts gain adequate immune protection, leading to extinction (stochastic fade out) [47–50], or rebounds if the viral population antigenically drifts sufficiently far from the ancestor strain to escape host immunity [Fig. 2(a)].

The escape probability depends on the total number of hosts, as we see by increasing the number of hosts in a one-deme case [Fig. 2(b)]. As $N_h \to \infty$, the viral survival probability in one deme approaches 1 because the virus has more time to develop crucial escape mutations before exhausting its supply of susceptible hosts.

Survival of the viral population also depends on the antigenic diversity of an outbreak, depicted schematically in Fig. 2(c). The variance of the density of infected hosts $n_1(x, t)$, $V_1(t) = \langle x^2 \rangle_1 - \langle x \rangle_1^2$, where averages are taken over variants in deme 1 according to measure $n_1(x, t)/N_1(t)$, gives an estimate of the diversity of variants in the population. Outbreaks with greater diversity have better access to novel mutations and are more likely to escape. Seeding outbreaks with Gaussian distributed initial $n(x, t = 0)$ with different values of $V(t = 0)$, we record the diversity at time $T$ of the outbreak peak, $V(T)$. We find a linear dependence between the escape probability and the diversity at the outbreak peak [Fig. 2(d)], showing that antigenic diversity influences viral survival.

#### B. Two demes

To see how the coupling of mobility and antigenic diversity influences viral escape, we turn to a two-deme system, with a constant migration rate $k$ between the two demes. The off-diagonal elements of the migration matrix $k \equiv K_{12} = K_{21}$ in Eq. (1) set the timescale over which the dynamics in the two demes equilibrate. We initialize our simulations with $N_0 = 100$ infected individuals in deme 1, as in the one-deme case, and none in deme 2.

At zero migration rate, individuals are unable to move between demes and the epidemiological dynamics in each deme are decoupled, replicating the case of a single deme discussed above. For infinite migration rates, spatial variation in the number of infected hosts is removed by rapid host movement and the combined population behaves as a single well-mixed population with twice the number of susceptible hosts. By calculating the proportion of trajectories in which the virus successfully escapes as a function of the spatial migration rate, we recover the expected behavior for small and large migration rates [Fig. 3(a)]. We also find that an intermediate migration rate maximizes the probability that the virus escapes host immunity, showing there is a preferred migration ratio that facilitates viral survival. The viral survival probability shows exactly the same dependence on the migration rate when we artificially forbid back
<figure><img src="image_4.png" /><figcaption>FIG. 2. For one deme, diverse pathogen populations better evade host immunity by accessing novel antigenic mutations. (a) Total number of infected hosts $N(t) = \int n(x, t)dx$ for 200 replicate realizations of the model dynamics with an isogenic initial antigenic distribution. All replicates follow a similar trajectory until the accumulation of population size noise leads to one of two outcomes: pathogen escape or extinction. The timing of the peak of the outbreak is denoted by $T$. (b) Escape probability as a function of the number of susceptible hosts. When the number of hosts is large, the probability of escape approaches 1. (c) Cartoon of the density of infected hosts in antigenic space $n(x, t)$ during the early outbreak. The area under the curve gives the total number of infected individuals $N(t)$ plotted in (a), while the variance of the antigenic density $V(t)$ gives a measure of the diversity of antigenic strains within the population. (d) Probability of antigenic escape as a function of antigenic diversity at the time of the outbreak peak. We generate trajectories with different levels of antigenic diversity by starting with a Gaussian infected density, with different initial variances. We measure the diversity at the peak of the outbreak $V(T)$ [time $T$ in (a)]. Larger initial variance leads to an increased diversity at the outbreak peak. We quantify this relationship with a linear fit $p = p_0 + m \times V(T)$ (dashed line) with slope $m$. The dependent variable in the linear regression $V(T)$ is measured from simulations and can be tuned by changing the diversity of the initial pathogen population $V(0)$. The other parameters are $\beta = 2.5$, $\alpha = 0$, $\gamma = 1$, $D = 0.01$, $\sigma = 2$, $M = 15$, $N_h = 2 \times 10^6$, and $r_0 = 3$. The analysis was performed over $10^4$ statistical replicates for (b) and $3 \times 10^3$ replicates were used for (d). The initial variances were linearly spaced between $10^{-2}$ and $10^{-1}$ with six initial variances being used in total.</figcaption></figure>

migration of infected individuals from deme 2 to 1, $K_{21} = k$ and $K_{12} = 0$ [Fig. 3(a), triangles], suggesting that the success of the secondary epidemics in deme 2 is key to the virus survival.

To better understand, quantitatively, the interplay of migration and antigenic escape, we consider the timescales of outbreak peaks in the two demes for the simpler case of symmetric migrations. Figure 3(b) shows that the total number of infected hosts in each deme shows a time delay $\Delta T = T_2 - T_1$ between their outbreak peaks. The distribution of infected hosts $n(x, t)$ in deme 2 at the outbreak peak is broader than in deme 1, signifying a higher diversity outbreak [inset of Fig. 3(b)]. The probability of antigenic escape for each migration rate is given by the probability of escape in each deme,

$p_1$ and $p_2$,

$$p_{\text{escape}} = 1 - (1 - p_1)(1 - p_2). \qquad (6)$$

In the short-time or low migration limit, we estimate the escape probability in deme 1, $p_1$, directly from one-deme simulations. The survival probability in deme 2, $p_2$, is the probability that the pathogen escapes due to its increased diversity, given that it spreads from deme 1 to deme 2. Due to the founder effect that gives an advantage to the very first migrating strains, the linear dependence of the escape probability on the antigenic diversity at the outbreak peak $V(T)$, as described in Fig. 2(d), should also hold for deme 2, so that we may write

$$p_2 = \{p_1 + m[V_2(T_2) - V_1(T_1)]\}p_{\text{spread}}, \qquad (7)$$
<figure><img src="image_5.png" /><figcaption>FIG. 3. Moderate host migration between two demes maximizes the probability of survival by amplifying the antigenic diversity of secondary outbreaks. (a) Probability of antigenic escape at different migration rates for symmetric 1 ↔ 2 (circles) and asymmetric 1 → 2 (triangles) migrations. The black dashed line is a prediction for the survival probability using the linear relation in Fig. 2(d) and Eq. (6). The horizontal lines show the results for the two limiting cases of zero migration rate (blue) and infinite migration rate (orange). The probability of survival at intermediate mobility attains a maximum near a rate of $k/γ ≈ 10^{-4}$. We do not report simulations of asymmetric migration for large $k$ because we neglect the time dependence of $N_h$, which should only be important for large $k$. (b) Number of infected hosts in each deme for symmetric migrations. The outbreak is seeded in deme 1 with 100 individuals and spreads through host migration to deme 2 at a time $ΔT$ in the future. The inset gives a comparison of the viral density of each outbreak at their respective peak $n_1(x, T_1)$ and $n_2(x, T_2)$. The outbreak in deme 2 has large variance, i.e., is more antigenically diverse. (c) Increased mean variance of the deme 2 outbreak compared to deme 1 at their respective peaks $⟨V_2(T_2) - V_1(T_1)⟩$ as a function of the average time difference $⟨ΔT⟩$ between the outbreak peaks. Both axes have been scaled by the cross reactivity $r_0^2$. The dashed line is the prediction $ΔV = 2DΔT$ [Eq. (8)]. For small migration rate, the time to seed an outbreak in deme 2, $ΔT$, is larger than for a large migration rate increasing the antigenic diversity, which increases the probability of antigenic escape as shown in Fig. 2(d). The blue solid line shows a prediction from the analytic expression in Eq. (F9), discussed in Appendix F. (d) Escape probability at different mutation and migration rates. The vertical axis shows the normalized escape probability, where the normalization factor for a given mutation rate is chosen so that the peak normalized escape probability is one. We see that the optimal migration rate is independent of the mutation rate even though the overall escape probability depends strongly on $D$. Only mutation rates with $10^{-3} < \text{escape probability} < 1$ are shown. The parameters are $β = 2.5$, $α = 0$, $γ = 1$, $D = 0.01$, $σ = 2$, $M = 15$, $N_h = 2×10^6$, and $r_0 = 3$. The migration rate in (b) is $k = 10^{-3}$. In (a), (c), and (d) analysis was performed over $10^4$ statistical replicates.</figcaption></figure>

where $m$ is the regression coefficient and $p_{spread}$ is the probability that the epidemic spreads from 1 to 2.

For each value of the migration rate, we record the average antigenic diversity at the outbreak peak in the two demes, $⟨V_1(T_1)⟩$ and $⟨V_2(T_2)⟩$, and the frequency with which outbreaks in deme 2 are seeded via migration $p_{spread}$. The expected overall escape probabilities calculated using Eq. (6) correctly predict the measured probabilities for all but very high mutation rates [black line in Fig. 3(a)], showing that the increased survival probability originates from an increased antigenic diversity in deme 2.

To understand why diversity in deme 2 is higher at its epidemic peak than in deme 1, we analyze how the variances $V_1(t)$ and $V_2(t)$ are predicted to evolve according to Eq. (1). During the first epidemic peak, immune pressure may be neglected and viral evolution in deme 1 is dominated by pure diffusion so that $V_1(t) ≈ 2Dt$. As soon as deme 2 receives some infected hosts from deme 1, its own diversity quickly tracks that of deme 1 when the migration is high enough, $V_2(t) ≈ V_1(t)$, so that

$$V_2(T_2) - V_1(T_1) \approx 2D\Delta T. \quad (8)$$
<figure><img src="image_6.png" /><figcaption>FIG. 4. Time dependence of the fixation index in the case of two demes. The solid lines show the average value over 200 replicate simulations and the shaded areas show the standard error.</figcaption></figure>

This prediction shows good agreement for small values of $\Delta T$, corresponding to large values of the migration rate $k$ [Fig. 3(c)]. Decreasing the migration rate increases the time delay $\Delta T$, increasing antigenic diversity and the probability of survival of the virus in deme 2. As the migration rate becomes very low and $\Delta T$ large, the secondary epidemic in deme 2 is seeded by only a few migrating hosts, leading to a collapse in the diversity at its peak. The effect may be accounted for qualitatively using simplifying approximations (see the Appendix A), yielding analytic or semianalytic predictions for the diversity at the epidemic peak [solid line in Fig. 3(c)]. Combining these predictions with Eq. (7) allows us to derive an expression for the optimal migration rate, which decreases as a function of the magnitude of the epidemic peak. The feedback between migration rate and escape probability is entirely mediated by the outbreak delay time $\Delta T$, which should be independent of the mutation rate $D$. In Fig. 3(d) we perform sweeps over both the mutation and migration rates and confirm that the optimal migration rate is insensitive to the choice of mutation rate. Outside of a small range around $10^{-2}$, the escape probability is either very close to 1 (high $D$) or very close to 0 (low $D$), making it difficult to see the effect of migration in our simulations. Within the range of interest, migration rate controls the sharpness of the transition, with the fastest rise occurring for the optimal migration rate.

To quantify the importance of the spatial structure of strain diversity, population geneticists often measure the fixation index of a fragmented population, which compares the diversity across the whole system to the average diversity in each deme. When computed as a function of time for a two-deme system, that fixation index reveals a period where the variances of strains in the two demes are markedly different (Fig. 4; see also Appendix D). This difference subsides at the end of the simulation as migration homogenizes the outbreaks. This confirms that the diversity boost created by spatial fragmentation is only transient. The magnitude of this effect depends nonmonotonically on the migration rate: There exists an intermediate migration rate for which it is highest. That rate is much smaller than the rate at which survival probability peaks, suggesting that the fixation index does not capture important aspects of immune escape.

Altogether, the analysis of the two-deme case shows that the optimal migration should be slow enough to produce a delay in the outbreak dynamics but fast enough so that the diversity of the initial epidemic in deme 1 can be inherited by deme 2. This effect entirely relies on the interplay of the migration and mutation rates.

### C. Linear networks with many demes

To investigate networks with more than two demes and the role of subsequent deme outbreaks, we first consider the simple topology of $N_D$ demes arranged with linear connectivity, each deme with $N_h$ hosts. We assume symmetric rates normalized by the number of outgoing connections a deme has such that the rate out of any deme sums to $k$. We observe a similar tradeoff as in the two-deme case [Fig. 5(a)], with an optimal migration rate achieving the balance between access to more hosts and the ability for additional demes to give extra chances of survival to the virus. Increasing the number of demes increases the escape probability at all migration rates, since it simply increases the number of hosts. This is clear in the limit of large migration rate ($k \to \infty$), which reduces to a single large deme of increasing size $N_D \times N_h$.

To disentangle the impact of deme structure from that of the number of hosts, we repeat the analysis but with constant total population size $N_{tot}$ and deme sizes that decrease accordingly with the number of demes as $N_h = N_{tot}/N_D$ (Fig. 5). Despite this normalization, which decreases the overall survival probability, adding demes still increases viral survival at intermediate migration values [Fig. 5(b)]. The optimal migration rate for viral survival shifts slightly towards larger values as the number of demes increases. This result is explained by the number of hosts per deme decreasing, which makes survival in individual demes harder and rescue through migration more important. Note that the large migration limit coincides with the single-deme result with $N_{tot}$ hosts, as expected.

The benefit of adding more demes even as the total population size is kept constant emphasizes the importance of inheritance of diversity between demes. Diversity builds faster within a deme when the viral population is pressured by the immune system. When a new deme is seeded with the virus, there is no established immune system and so not much diversity is generated inside it. Instead, diversity is imported from the previous deme, where it is already large, as additional infected hosts arrive.

For fixed numbers of hosts per deme [Fig. 5(a)], the optimal migration rate for viral survival is the same for all network sizes. As we argue in the Appendix F, the location of the peak decreases with the maximal number of hosts infected during the initial outbreak, $N_{max}$. That peak is a property of each newly infected deme. It depends not on the overall number of demes but only on the number of hosts in the infected deme. When we add demes of constant size $N_h$ to the system, $N_{max}$ remains similar for each deme and so does the optimal migration rate. By contrast, when we increase the number of demes while keeping the total population constant, we expect the migration rate to shift to higher values as $N_h$, and thus $N_{max}$, decreases, as observed in Fig. 5(b).

In Fig. 5(c) the importance of viral diversity for survivability in a larger deme system is revealed. The variance of the
<figure><img src="image_7.png" /><figcaption>FIG. 5. Increase of survival probability with the number of demes. The survival probability is plotted as a function of the migration rate for networks with $N_D$ demes connected in a linear topology, with (a) a fixed number $N_h = 10^6$ of hosts per deme and (b) a fixed total population $N_{tot} = 10^6$ and $N_h = N_{tot}/N_D$ hosts per deme. In both cases, symmetric migration rates are normalized by the number of outgoing connections a deme has such that the rate out of any deme sums to $k$. (c) Mean variance of the viral antigenic density within each deme when the number of infected people in the deme is maximized in a five-deme system with line topology; the same system producing the red line in (a) was used. Different demes along the line are indicated by different symbols and color shades.</figcaption></figure>

viral density is measured for each deme in a five-deme line topology system. The last deme achieves the largest variance in the system, demonstrating the role of adding more demes in increasing viral diversity and therefore survival. Diversity in each deme is maximized at the migration rate at which survival probability also reaches a peak. This optimum occurs when the variance inherited by each successive deme is the same and the total variance within the system is also maximal. Consequently, the total diversity is optimized at a migration rate which allows for the same inherited diversity in all demes, leading to the best chance of survival.

Overall, we conclude that separating hosts into many demes helps the viral population survive.

### D. Influence of network topology on viral escape probability

To explore the impact of topology beyond linear networks, we consider all unique four-deme network topologies [Fig. 6(a)], where uniqueness is defined with respect to network structure and the identity of the initial outbreak deme. We assume symmetric migration rates normalized by the number of outgoing connections such that the overall migration rate out of each deme is $k$. Figure 6(b) shows the survival probability as a function of migration rate for these topologies. As in the previous cases, there exists an optimal migration rate at which the probability of survival is maximal, and the low and high migration rate limits reduce to the single-deme case with $N_h$ and $N_{tot} = 4N_h$ hosts. The optimal migration rate

<figure><img src="image_8.png" /><figcaption>FIG. 6. Effect of network topology on virus survival probability. (a) All four-deme networks. White triangles indicate the deme/node with the initial outbreak. The color represents the value of the centrality of the outbreak node. Open symbols in the lower right corner of each colored box refer to the markers used in (b), specifying the network topology used for the simulation of each curve. Symmetric migration rates are normalized so that the total outgoing migration rate of each deme is always $k$. (b) Viral survival probability as a function of migration rate for all possible four-deme networks. The inset shows the peak height as a function of centrality. The simulation parameters are $\beta = 2.5$, $\alpha = 0$, $\gamma = 1$, $D = 0.01$, $r_0 = 3$, $M = 15$ $\sigma = 2$, and $N_h = 10^6$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Migration rates for real networks given flight travel data. (a) The top 20 most frequently traveled airports from the data are displayed in a network graph with color gradient representing the closeness centrality of each deme. (b) Migration rates are calculated for different real networks and are put into context with the four-deme line topology migration rate which maximizes the survival probability. How these rates were obtained is outlined in Appendix G. Data for the number of flights within designated regions are collected from [52]. (c) An inverse relationship between closeness centrality and probability of survival of the virus exists within the real network data. The parameters for simulation are $\beta = 2.5$, $\alpha = 0$, $\gamma = 1$, $D = 0.0042$, $r_0 = 3$, and $M = 15$ $\sigma = 2$. Host population sizes in each deme reflect the metropolitan area population of the city where each airport resides [53].</figcaption></figure>

for viral survival is relatively insensitive to network structure, since it is mostly driven by the deme size, as discussed earlier. However, the value at the peak depends on the particular network topology.

To identify what features of the topology drive this dependence, we consider closeness centrality, a quantity that was proposed to quantify the importance of a node $x$ in a network and defined as the inverse average distance to all other nodes in the system [35]

$$C(x) = \frac{N - 1}{\sum_{y=1}^{N-1} d(y, x)}, \qquad (9)$$

where $N$ is the total number of nodes in the network and $d(y, x)$ is the length of the shortest path between nodes $x$ and $y$. We colored networks in Fig. 6(b) and the survival probability curves in Fig. 6(a) according to the closeness centrality of the outbreak node. The smallest closeness centrality value, represented by the line topology (in green), results in the highest maximum survival probability. Networks with nodes with high closeness centralities have lower survival probability peaks [inset of Fig. 6(a)], with the fully connected network having the lowest. Outbreaks with the same centrality have very similar curves.

Two aspects contribute to this effect: the network topology itself and the location of the outbreak. Lower connectivities enhance virus survival, as do outbreak demes that are isolated. Network topology introduces an exploration-exploitation tradeoff to viral survival. Outbreaks originating in

nodes that are well connected gain access to more demes, and therefore susceptible hosts, early on. However, as the outbreak spreads quickly, only the diversity developed in the first deme is transmitted to other demes. Sequential discovery of new hosts in less connected networks allows the viral population to accumulate diversity via sequential range expansions. In the two-deme case, access to new populations at a delayed time drives an increase in viral diversity. In larger networks, the structure of the network modulates the number of times the virus can exploit diversity accumulation due to this time delay. Outbreaks in well-connected nodes expand to new populations simultaneously, while the linear structure of $N_D$ demes allows $N_D - 1$ temporally ordered range expansions that optimally exploit each deme.

### E. Where do real travel networks lie?

To determine if the relationship between closeness centrality and viral survival probability can be observed in real network data, we analyze an airplane travel data set [51]. The full data set is composed of 3632 airports and the number of people who travel between airport pairs in 2011. We select the top 20 most frequently traveled to airports and plot the network connectivity graph with color coding according to the closeness centrality found for each node [Fig. 7(a)]. The closeness centrality is calculated first by designating local distances between any two cities as the inverse of the number of people traveling between them. Using these distances as the directed weights in the network, Dijkstra's algorithm is used
<figure><img src="image_10.png" /><figcaption>FIG. 8. Analysis of migration rate heterogeneity. The relationship between viral survival probability and the closeness centrality is measured for (a) the top 30 airports experiencing the most travelers arriving and departing, two bootstrap cases (b) and (c) where migration rates in the adjacency matrix are randomly selected from the distribution of migration rates from the top 30 airports, and (d) bootstrap simulation with half the migration rates drawn from the same distribution but symmetry in the adjacency matrix is enforced. The parameters are the same as in Fig. 7, with $N_h = 10^7$ for all demes in (b)–(d).</figcaption></figure>

to find the shortest path between the origin deme and all other demes in the system. More explicitly, from Eq. (7), $d(x, y) = \min(\sum_{zz'} \text{edges in path from } x \text{ to } y \ell_{zz'})$, where $\ell_{zz'}$ for each edge is the inverse of the number of people traveling along that edge. This places more strongly coupled airports closer together in network space. The network structure shows expected trends: The four Chinese airports (PEK, PVG, CAN, and CTU) are strongly connected, as are Mumbai (BOM) and Delhi (DEL) and the three U.S. airports (ATL, LAX, and DEN). We interpret each airport node of the graph as a deme.

To understand how our model predictions connect to real-world networks, we estimate travel rates based on the 20-deme network in Fig. 7(a), as well as other data sets [52]. Our estimates include only air travel and are likely lower than the real values which include other modes of transportation.

We then compare these rates to the optimal migration rates obtained by the model for the four-deme line topology, explored previously [Fig. 7(b)]. While topologies of real networks are far bigger and more complex than our simulations allowed us to explore, our previous results suggest that the location of the maximum depends only weakly on the details of the network topology or size. The comparison between the predicted optimal rate and the empirical estimates gives good agreement, suggesting that current travel volumes maximally improve a virus's chances of survival.

To delve more finely into the network structure, we use the adjacency matrix produced by the 20 airport system outlined in Fig. 7(a) to simulate and calculate the viral survival probability for outbreaks initiated at each of the 20 airports. The migration rates used reflect the number of people moving between airports and are normalized by the metropolitan area population size of the origin deme [53]. The viral escape probabilities are then plotted as a function of the closeness centrality of the airport where the outbreak originates. The viral survival probabilities are negatively correlated with the closeness centrality [Fig. 7(c)], meaning the more remote the outbreak deme, the greater the probability of survival. The same negative relationship between probability of survival and closeness centrality of the outbreak deme holds if we extend the analysis to the top 30 most trafficked airports [Fig. 8(a)]. Therefore, we observe the same relationship in real-world networks as in our smaller toy networks. The

network graph and resulting ranking of demes by closeness centrality predict that an outbreak originating in the deme with the lowest closeness centrality, the only South American airport São Paulo (GRU) in our 20 top airport reduced data set, has the highest probability of survival. Conversely, a virus originating in Hong Kong (HKG), which is well connected to other destinations around the world, benefits less from the ability to migrate.

To investigate how the distribution of migration rates and the distribution of population sizes influence the relationship between survival probability of the virus and the closeness centrality of the origin deme where the outbreak is initiated, we simulate two data sets of 30 airports with two different randomly sampled adjacency matrices from migration rate distributions obtained from real data and all initial population sizes for each of the 30 demes set to the same $N_h$ value (see Appendix 1 for details). Results are shown in Figs. 8(b) and 8(c). The negative trend is absent in those two synthetic networks, revealing that some structure is required for closeness centrality to act as an indicator of viral survival [Figs. 8(b) and 8(c)]. Unlike those synthetic networks, which sample all the rates independently, the real airport network is roughly symmetric. To test the toll of symmetry of the adjacency matrix, we randomly sampled migration rates as above but imposed matrix symmetry (see Appendix I). This procedure restores the negative correlation between closeness centrality and survival probability [Fig. 8(d)], indicating that what is important is that the fluxes in and out of each deme are roughly equal.

Finally, to understand the impact of population diversity, we designed another synthetic network where all migration rates were set to $10^{-4}$ and population sizes were sampled from the distribution of all metropolitan areas in the data. This analysis, shown in Fig. 9, shows no relationship between survival probability and population size. The flux out of a given deme in this simulation is governed by its population size, but since the adjacency matrix codes for a fully connected network, flux into each deme is determined by all other deme's populations. This means that a large deme in the system will most likely have a larger migration rate out than in. This imbalance erases the relationship between survival probability and deme size, similarly to what we observed for closeness centrality [Figs. 8(b) and 8(c)].
<figure><img src="image_11.png" /><figcaption>FIG. 9. Analysis of population size heterogeneity. Survival probability is plotted against population size. A negative trend reveals that outbreaks starting at demes with larger flux (larger population size) have a worse chance of survival. However, no relationship is identified.</figcaption></figure>

#### IV. DISCUSSION

When examining the interplay between viral evolution and host spatial migration we found that intermediate migration rates maximize viral diversity, favoring antigenic escape. This result stems from coupling a stochastic property of evolution, i.e., diversification, and exploiting susceptible individuals due to migration. Large migration rates allow diversity to be transmitted between infected demes, as newly infected populations exert only weak or absent selection and multiple migration events allow multiple strains to simultaneously grow. This coupling leads to greater diversity and a higher probability of viral survival and ultimately antigenic escape. This process is an example of an exploration-exploitation tradeoff: Low migration rates hamper the virus's search for new pools of susceptible individuals whereas high migration rates result in infection with similar strains, decreasing diversity and future antigenic escape. This tradeoff depends on the ratio of timescales for host immunity and migration. Low migration rates that result in exploiting a deme for long times allow the host immune systems to gain cohort immunity and eliminate viral diversity through selection. Long-term survival of the virus therefore exploits short-term coupling of host-induced mutations and migrating at peak diversity. This process requires dwell times in a population that allow for the strain diversification but avoid low migration rates that would lead to a decrease in diversity before migration.

The existence of an optimal migration rate for pathogen survival has been observed in many models and attributed to different effects. In our model, the optimal rate exploits the diversity buildup due to host-driven viral mutations, which are not included in models of long-term and large-scale outbreaks [33,54]. Matching allele models, originally introduced to study plant pathogen evolution, incorporate host/parasite migration, evolution, and demographics via a multicompartmentsal Lotka-Volterra style approach [55]. However, they assume fixed diversity. There are *n* possible alleles for a single locus in both parasite and hosts, in which the host is immune to the “matching” parasite strain but susceptible to all others. Within this class of models, allowing only for host migration, Gandon

*et al.* [55] reported an intermediate optimal value of host migration rates for host resistance. High pathogen persistence probabilities at intermediate pathogen migration rates have also been observed in two-compartment susceptible-infected-susceptible (SIS) models without an evolutionary component [33,56,57]. Unlike these results, the optimal migration rate for viral survival in our model appears also in small networks, whereas in the work of Hisi *et al.* [57] the optimal migration rate disappears in networks of ten or less demes.

The migration of susceptible hosts on the network has been discussed in terms of favoring one of two alternative rescue effects: evolutionary [55] or demographic [57]. The evolutionary rescue effect relies on gene flow, where the genetics of a population is affected by the genetic makeup of neighboring populations. This allows migrating strains that are poorly adapted to a given host environment to survive by migrating to new niches. It also allows well-adapted populations to spread their genes to other regions, accelerating recovery after a shock [41]. This process is distinct from the demographic rescue effect, which creates short-lived refuges for viruses by randomly bringing a large number of susceptibles together for a period of time. The viral population can survive in this ephemeral host population until it reinfects other regions that have by then lost their immunity. Our observations rely on evolutionary rather than demographic rescue mechanisms since removing back migration does not affect viral survival probabilities for a wide range of migration rates and viral survival correlates with large strain diversity. In general, they reinforce the importance of evolutionary dynamics in metapopulation persistence, in line with [40,41,58].

A number of analytic studies have focused on ecological models without evolution, describing demographic rescue effects in two-deme SIS systems. Meakin and Keeling [59] and Keeling and Rohani [60] use a moment closure approach to derive the correlation functions of the outbreak sizes in the demes. These steady-state results find highly correlated outbreaks (and therefore large maximum outbreak sizes) at intermediate coupling strength. (Both models use a different type of coupling: Rather than permanently migrating, individuals “commute” to the other deme for a fixed amount of time and then return home.) Of course, survival probability and population size are not the same thing: Our outbreaks peak after infecting most of the population, and then surviving simulations tend to recover and level off at about 10%-20% $N_h$ in each deme, roughly independent of migration rate. We also find that the outbreak peaks occur at the same time in different demes at high migration rates [see Fig. 3(c), where $\Delta T$ goes to 0 at high migration rates], with a peak asynchrony (large $\Delta T$) at intermediate times. However, the composition of the outbreaks in each deme is not the same, which allows us to draw evolutionary conclusions from this observation. The outbreak in deme 2 is much more diverse than in deme 1 (see Fig. 10 in Appendix C for histograms of variance for each deme). Since diversity drives antigenic escape [Fig. 3(b)], comparison with this class of models also leads us to conclude that evolution is responsible for the optimal migration rate for viral survival. At intermediate migration rates, the escaping strain appears in the more diverse deme 2, while at low migration rates it appears in deme 1, where the outbreak began (see Fig. 11 in Appendix C).
Evolutionary models without migration have predicted a nonlinear increase in the speed of antigenic evolution as a function of the fraction of immunocompromised hosts in the population [37]. This result is due to hosts that recover more slowly, infecting more people, as they have more time to transmit the infection. This causes an effective increase in the growth rate of the virus, increasing the probability of acquiring an antigenic mutation. Immunocompromised hosts play a role similar to migration in our model: When a new deme is infected for the first time, the virus has access to a large number of hosts with susceptible immune systems, increasing the growth rate of the virus, but only in a transient way. While we do not systematically study the effect of base virulence $\beta$ in our model, previous theoretical work has shown that shorter, more virulent outbreaks tend to produce less genetic diversity than longer, less infectious ones [36]. In other words, the virus needs time to exploit mutations, and time is in short supply in fast, sudden, and one-time outbreaks that occur with high $\beta$. Comparison to this class of models allows us to draw conclusions about the evolutionary impact of population structure in space.

The topology of the network determines the probability of survival of the viral population. Outbreaks in demes with low closeness centrality are more likely to escape the host immune systems than ones starting in high closeness centrality demes. This means that the way a virus disperses has evolutionary consequences, at least when evolution happens quickly. Although we focused on small networks, these effects probably hold in large networks. If this is the case, not all initial conditions of the simulation are equal. In the case of randomly chosen initial conditions over many replicates, not all replicates contribute equally to observed mean probabilities. Depending on how extreme this effect is, some initial outbreak locations may dominate viral survival probabilities on large networks. In the context of protecting crops from dispersing pests, it has been noted that the degree of the outbreak node is correlated with the size of the outbreak and the speed of spread [61,62].

The role of topology for viral survival also suggests that multiple pairwise mobility restrictions between cities or airports can synergistically strengthen and result in a stronger global effect. This has been reported when implementing spectral control strategies [63,64], which selectively modify traffic between specific nodes to control epidemics. In the long-time limit, SIS models on networks have shown that the maximal eigenvalue of the network determines the exponent of the rate with which a particular virus goes extinct [35,64]. Our model complements these results with analysis of the transient dynamics. The networks that favor antigenic escape in the short term are not necessarily the same as those that favor long-term survival. For example, the linear network topology guarantees the largest probability of short-term antigenic escape, but has a small maximal eigenvalue, which makes viral survival more difficult at long times. The network that makes viral survival most difficult at short times, the fully connected network, has a large maximal eigenvalue and offers the best conditions for long-term viral survival. Therefore, the traffic restriction strategies which are most effective will be different for emergent and endemic pathogens.

In the case of emergent pathogens, the applicability of our results will depend on the scale considered. Figure 7 shows that depending on which set of airport data one considers, the effective migration rate is not the same. All reported migration rates result in viral survival probabilities that are close to the maximum [Fig. 7(b)]. This result could imply that reducing viral survival could be achieved by either strongly increasing or decreasing mobility. However, the value of viral survival also depends on the number of hosts [Fig. 5(c)]. If the number of hosts is large, increasing migration rates does not decrease the viral survival probability. Additionally, the width of the peak in terms of migration rates means obtaining a significant effect would require changing mobility by several orders of magnitude, which is no small task for policymakers. For example, Ref. [65] report that cell phone mobility data show that the “lockdown” policy caused a 54% reduction in mobility in the Chicago metro area in the first week of April 2020. Additionally, while decreasing mobility is achievable, increasing mobility is harder to implement. Finally, driving a virus to extinction is not the primary goal of policy. It is more important to reduce the total number of cases and to prevent an outbreak large enough to overwhelm the healthcare system.

Eco-evolutionary feedbacks remain critical to understanding long-term epidemiology. We have examined a density-dependent selection effect (very common strains within a deme experience the most immune pressure, neglecting the lag time for collective immunity to catch up), but global frequency-dependent effects can be important as well. It also remains to be seen how short-time selection biases towards diversity impact mid- and long-term phylodynamics. Such a question requires speciation dynamics and therefore a higher-dimensional antigenic space.

Low-dimensional antigenic spaces are commonly used by experimentalists and theorists alike [18,66-68]. While it has been suggested that this low dimensionality could be an artifact of random walks [69], canalization of evolutionary trajectories in the long-time endemic state suggests that the phenomenon is not an artifact, but rather is caused by the propagation of an antigenic wave [19]. However, even in a canalized state, trajectories are not predictable beyond a persistence time that decreases with dimension of the antigenic space [23].

The effective antigenic dimension of emerging pathogen coevolution has been much less discussed. Our results may be sensitive to the choice of dimension, as higher dimensions should generically increase antigenic diversity and increase the probability of antigenic escape. However, we expect the qualitative structure of the peak at intermediate migration rates to hold for any dimensionality, because our analytic approach (see Appendix A) readily extends to any dimension in the Gaussian approximation. We therefore expect a higher, but still nonmonotonic escape probability as a function of the migration rate.

## ACKNOWLEDGMENTS

This work was carried out as a summer school project during the Les Houches summer school on Theoretical Biology. A.M.W. and T.M. thank Alain Barrat for useful discussions. This work was supported by the European Re-
search Council consolidator Grant No. 724208 (A.M.W., T.M.), the Agence Nationale de la Recherche Grant No. ANR-19-CE45-0018 “RESP-REP” (A.M.W., T.M.) and the Chan Zuckerberg Theory Initiative (A.M.W.). AI-R received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 847718. D.W.S. acknowledges the MIT Engaging cluster for providing computational resources and support. N.B. acknowledges funding from the EMERGENCE(S) grant program of the city of Paris. This publication reflects only the authors’ views. The relevant funders are not responsible for any use that may be made of the information it contains. C.B. thanks funding from the NSF GRFP.

## DATA AVAILABILITY

All codes used to perform simulations and generate figures can be found in [70].

## APPENDIX A: OLIGOMORPHIC DYNAMICS

We adapt the oligomorphic dynamics (OMD) approximation [71,72] to study the size and antigenic variance of the initial outbreak. For simplicity, we neglect noise and work in two demes. The dynamics of the densities in deme i are given by Eqs. (1)–(4), where we define the fitness as $F_i(x, t) = \beta S_i(x, t) - \alpha - \gamma$ and, for simplicity, we neglect noise and work in two demes. We assume that there are $N_0$ individuals infected at time $t = 0$, all with strain $x = 0$ in deme 1, which gives the initial condition $n_1(x, t = 0) = N_0\delta(x)$ and $n_2(x, t = 0) = 0$. We obtain an equation of motion for the total number of infected individuals $N_i(t) = \int dx n_i(x, t)$ by integrating Eq. (1):

$$ \frac{dN_i}{dt} = \langle F_i \rangle N_i + k(N_j - N_i). \quad (A1) $$

The total infected population grows with the antigen-averaged fitness and migration causes the infected numbers in different demes to equilibrate on a timescale of $k^{-1}$.

We define the densities of viral variants and immune protections as

$$ \phi_i(x, t) = \frac{n_i(x, t)}{N_i(t)} \quad (A2) $$

$$ \psi_i(x, t) = \frac{h_i(x, t)}{H_i(t)}, \quad (A3) $$

with $H_i(t) = \int dx h_i(x, t)$. Applying Eqs. (A2) and (A3) to Eq. (1), we obtain

$$ \begin{aligned} \frac{\partial \phi_i}{\partial t} &= \frac{1}{N_i} \frac{\partial n_i}{\partial t} - \frac{n_i}{N_i^2} \frac{dN_i}{dt} \\ &= (F_i - \langle F_i \rangle) \phi_i + D \frac{\partial^2 \phi_i}{\partial x^2} + \frac{k N_j}{N_i} (\phi_j - \phi_i). \end{aligned} \quad (A4) $$

The first term tells us that regions of antigen space with fitness larger than average will increase in frequency, while regions with fitness less than average will decrease in frequency. The migration term has an $N_j/N_i$ prefactor. If deme j has a much larger infected population than deme i, the antigenic distribution in deme i will rapidly relax to the antigenic distribution

of deme j, inheriting all the genetic diversity from deme j almost immediately. This rapid equilibration of the antigenic distributions among demes explains the enriched diversity vs outbreak size relationships for secondary outbreaks we see in simulation.

We now assume $\phi(x, t)$ is Gaussian and can thus be characterized completely by its mean $m_i(t) = \langle x \rangle_i$ and variance $V_i(t) = \langle (x - \langle x \rangle_i)^2 \rangle_i$. We assume the distribution of $x$ is sharply peaked about the mean value $m_i(t) = \langle x \rangle_i$ and expand the fitness around the mean

$$ \begin{aligned} \langle F_i(x) \rangle_i &\approx \langle F_i(m_i) + F_i'(m_i)(x - m_i) \rangle_i \\ &\quad + \frac{1}{2} F''(m_i)(x - m_i)^2 + \dots \rangle_i \\ &= F_i(m_i) + \frac{F_i''(m_i)}{2} V_i^2 + \dots, \end{aligned} \quad (A5) $$

where we neglect higher-order terms than the variance $V_i(t) = \langle (x - \langle x \rangle_i)^2 \rangle_i$. This approximation is valid since for typical simulation parameters $V/r_0^2 \sim 2Dt/r_0^2 \sim 10^{-2}$, which renders even the higher-order terms small compared to first- and second-order terms. The OMD approximation for the total infected number gives

$$ \frac{dN_i}{dt} = \left( F_i(m_i) + \frac{F_i''(m_i)}{2} V_i \right) N_i + k(N_j - N_i). \quad (A6) $$

To close the equation we need the dynamics for the mean and variance $m_i$ and $V_i$, respectively; $m_i = 0$ at all times based on symmetry arguments. Multiplying Eq. (A5) by $x^2$ and integrating by parts, we find

$$ \begin{aligned} \frac{dV_i}{dt} &= \int x^2 \left( (F_i - \langle F_i \rangle) \phi_i + D \frac{\partial^2 \phi_i}{\partial x^2} + \frac{k N_j}{N_i} (\phi_j - \phi_i) \right) dx \\ &= \mathrm{Cov}(x^2, F_i(x)) + 2D + k \frac{N_j}{N_i} (V_j - V_i). \end{aligned} \quad (A7) $$

This is a generalization of the Price equation to include diffusive spreading due to mutations and the rapid distribution equilibration due to migration.

## APPENDIX B: DERIVING $\Delta V = 2D\Delta T$

For the case of two demes, the covariance term in Eq. (A7) scales with $V_i^2$ at early times when $V_i$ is still small and can thus be neglected, giving equations for the variance in the OMD approximation

$$ \frac{dV_1}{dt} = 2D + k \frac{N_2}{N_1} (V_2 - V_1), \quad (B1) $$

$$ \frac{dV_2}{dt} = 2D + k \frac{N_1}{N_2} (V_1 - V_2). \quad (B2) $$

If the outbreak starts in deme 1 and spreads to deme 2 at $\Delta T$, at small times we have $V_1 \approx 2Dt$ and the solution of

$$ \frac{dV_2}{dt} = 2D + k \frac{N_1}{N_2} (2Dt - V_2), \quad (B3) $$

$$ V_2(\Delta T) = 0 \quad (B4) $$

is

$$ V_2(t) = 2D \left[ t - \Delta T \exp \left( -k \int_{\Delta T}^{t} \frac{N_1(t)}{N_2(t)} dt \right) \right]. \quad (B5) $$
<figure><img src="image_12.png" /><figcaption>FIG. 10. Histograms of diversity at outbreak peak in each deme show that at optimal migration rates, deme 2 is markedly more diverse than deme 1. The histograms are count normalized, so low migration rates which rarely seed secondary outbreaks produce very few samples of diversity in deme 2.</figcaption></figure>

If the infected number of individuals in deme 1, $N_1$, peaks at time $T_1$ and the infected number in deme 2 peaks at time $T_2 = T_1 + \Delta T$, the antigenic variance at the outbreak peak in deme 2 is

$$V_2(T_1 + \Delta T) = 2DT_1 + 2D\Delta T \times \left[ 1 - \exp \left( -k \int_{\Delta T}^{T_1+\Delta T} \frac{N_1(t)}{N_2(t)} dt \right) \right]. \quad (B6)$$

For large $k$, the spread to deme 2 happens quickly, $\Delta T \approx 0$, such that

$$V_2(T_2) \approx 2D(T_1 + \Delta T) = V_1(T_1) + 2D\Delta T. \quad (B7)$$

## APPENDIX C: ANTIGENIC DIVERSITY IN EACH DEME

We have argued in the main text that antigenic diversity drives our results. Here we include more detailed simulation results from the two-deme case to demonstrate that diversity is systematically higher in the second deme. We compare diversity between demes at comparable points in the development of the outbreak in each deme, i.e., the diversity of each deme at the peak of the outbreak in each deme. This amounts to comparing the within deme diversity with the time lag of $\Delta T$.

Figure 10 shows the histogram of antigenic diversity in each deme at outbreak for a sweep of migration rate $k$. At low migration rates (bottom) there are very few outbreaks in deme 2, so the deme 2 histogram is very sparse. As migration rate increases and becomes close to the optimal migration rate for viral survival, we see that deme 2 tends to have more diverse outbreaks, with the mean of the distribution shifting right as well as its width growing to close to fourfold compared to deme 1. Deme 2 is much more likely to see very diverse outbreaks than deme 1. Increasing $k$ further, rapid mixing between the demes takes over and the distributions become

<figure><img src="image_13.png" /><figcaption>FIG. 11. Probability of escape occurring in each deme for various migration rates in the two-deme case. The size of the red region gives the probability that the antigenically escaped strain first appears in deme 2 as opposed to deme 1. We define an antigenically escaped strain as one that is sufficiently far from the ancestor strain at $x = 0$. In deme $i$, the time of emergence is the first time $t$ such that $\int_{|x|>x^*} n_i(x, t)dx > 100$, with $x^*$ (the critical antigenic distance) computed as in Ref. [38]. If the new strain emerges in both demes within two time units of one another, we say they appear at the same time. We highlight that at intermediate migration rates deme 2 is the primary source of escape strains due to its amplified antigenic diversity. The parameters are the same as in Fig. 3.</figcaption></figure>

more similar. Here invasion is too rapid to exploit deme 1 to gain diversity and the virus immediately explores both demes.

The correlation of higher diversity with antigenic escape is demonstrated by recording which deme produces an escape variant first. We use OMD arguments to set an escape distance $x^*$ such that immunity profile created in response to an outbreak centered at 0 has largely decayed. We define the escape time $T_e^i$ for each deme $i$ as the first time a virus appears in that region of antigenic space. We then compare the escape times in each deme to see, over many runs, where the escape is most likely to occur (Fig. 11). At low migration rates, deme 1 has escape events more often than deme 2, as deme 2 is more rarely populated. At intermediate migration rates, the more diverse deme 2 has escape events more often than deme 1, and at high migration rates the two demes show similar numbers of escape events.

## APPENDIX D: FIXATION INDEX

A classic way to assess the importance of spatial structure on strain diversity is to compute the fixation index $F_{ST}$, which compares the proportion of genetic diversity due to local and global variability

$$F_{ST}(t) = \frac{V_{within}(t)}{V_{total}(t)}, \quad (D1)$$

where $V_{within}$ is the antigenic variance within a typical deme and $V_{total}$ is the antigenic variance across all demes in the network. The number of infected individuals is different between the two demes, so we compute $V_{within}(t)$ as a weighted average

$$V_{within}(t) = \frac{N_1(t)V_1(t) + N_2V_2(t)}{N_1(t) + N_2(t)}. \quad (D2)$$
The total antigenic variance $V_{\text{total}}$ is found by combining the antigen-space densities across all demes

$$n_{\text{total}}(x, t) = \sum_{i} n_{i}(x, t). \qquad (\text{D3})$$

The combined viral density is then used to compute the total variance of the antigenicity $V_{\text{total}}$.

## APPENDIX E: NUMERICAL PROCEDURES

Our simulations produce solutions to the system of stochastic integro–partial differential equations (PDEs) defined by Eqs. (1)–(4), with all codes to produce simulation results and figures available from [70]. Our simulation method draws on previous work, mainly Refs. [73–75], which we review here. The key idea behind our numerical solver is the operator splitting method [76,77], where the terms in Eq. (1) are solved one at a time, sequentially for each time step. We split Eq. (1) into a deterministic and a stochastic piece and solve each individually. We solve our system on a uniform grid of antigenic points with grid spacing $\Delta x$ and periodic boundary conditions. The length of the antigenic space is $L$, so the antigenicity can take values $x = -L/2, -L/2 + \Delta x, \dots, L/2 - \Delta x$. We assume periodic boundary conditions so that the antigenic site $x = L/2$ is identical to $x = -L/2$, but other choices of boundary condition (such as no flux) are possible. We choose values of $L$ large enough so that the infected density never interacts with the boundary of the antigenic space. As in the continuum model, we will call the infected density in deme $i$ at antigenic coordinate $x$ and time $t$, $n_i(x, t)$, the immune density $h_i(x, t)$, the cross-reactive immune protection $c_i(x, t)$, and the susceptibility $S_i(x, t)$.

Time is likewise discretized by steps of length $\Delta t$; for the deterministic step of the algorithm we use the explicit forward-time Euler method to update all densities. In a given time step $\Delta t$, we first perform a diffusion update to capture mutational changes. At each antigenic point $x$ we compute

$$\tilde{n}_i(x, t) = n_i(x, t) + \frac{D\Delta t}{\Delta x^2} [n_i(x + \Delta x, t) + n_i(x - \Delta x, t) - 2n_i(x, t)]. \qquad (\text{E1})$$

We then compute the cross-reactive immune protection $c_i(x, t)$ at each grid point $x$,

$$c_i(x, t) = \sum_{y} \exp \left( \frac{\min(|x - y|, L - |x - y|)}{r_0} \right) h_i(y, t) \Delta x, \qquad (\text{E2})$$

where the sum runs over $y = -L/2, -L/2 + \Delta x, \dots, L/2 - \Delta x$ and the minimum function accounts for the periodic boundary conditions on antigenic space. The cross-reactive protection $c_i(x, t)$ is then used to compute the susceptibility $S_i(x, t) = [1 - c_i(x, t)]^M$. We next update the immune and infected densities using the postmutation infected density $\tilde{n}_i(x, t)$,

$$h_i(x, t + \Delta t) = h_i(x, t) + \frac{\Delta t}{MN_h} [\tilde{n}_i(x, t) - \tilde{N}_i(t)h_i(x, t)], \qquad (\text{E3})$$

$$\tilde{n}_i(x, t) = \tilde{n}_i(x, t)[1 + \Delta t F_i(x, t)], \qquad (\text{E4})$$

where

$$\tilde{N}_i(x, t) = \sum_x \tilde{n}_i(x, t) \Delta x \qquad (\text{E5})$$

and

$$F_i(x, t) = \beta S_i(x, t) - \alpha - \gamma. \qquad (\text{E6})$$

We now turn to treat the stochastic piece of the dynamics. The characteristic feature of demographic noise is fluctuations that scale like the square root of the population size. Between any two grid points in antigenic space $x \to x + \Delta x$, we expect there to be $n_i(x, t)\Delta x$ infected individuals in deme $i$. We then expect that the noise amplitude should be proportional to $\sqrt{n_i(x, t)\Delta x}$. The factor of $\Delta x$ here is important; previous work [73–75] assumes a unit lattice spacing, which is appropriate for a stochastic PDE. Our model by contrast requires direct integration over the antigenic coordinate in order to compute immune response. This integration requires definition of an antigenic step size $\Delta x$, which we eventually make small. We begin by writing the stochastic differential equation (SDE) for the number of infected individuals between $x$ and $x + \Delta x$,

$$\frac{d}{dt}(n_i(x, t)\Delta x) = \sigma \sqrt{n_i(x, t)\Delta x} \xi_i(x, t). \qquad (\text{E7})$$

The process $\xi_i(x, t)$ is a white noise with correlation function

$$\langle \xi_x(x, t) \xi_j(y, s) \rangle = \delta_{ij} \delta_{xy} \delta(t - s). \qquad (\text{E8})$$

The SDE for the density at grid point $x$ is then

$$\frac{dn_i}{dt} = \frac{\sigma}{\sqrt{\Delta x}} \sqrt{n_i} \xi_i(x, t). \qquad (\text{E9})$$

The solution to this SDE is known [75]. We directly integrate from $t$ to $t + \Delta t$, where the initial condition is the result of the deterministic step $\tilde{n}_i(x, t)$. The distribution for $n_i(x, t + \Delta t)$ is

$$n_i(x, t + \Delta t) \sim \text{Gamma} \left( \text{Poisson} \left( \frac{2\tilde{n}_i(x, t)\Delta x}{\sigma^2 \Delta t} \right) \right) \frac{\sigma^2 \Delta t}{2\Delta x}. \qquad (\text{E10})$$

The argument of the Gamma distribution above is the scale parameter; the rate parameter is 1. Other methods of sampling the noise are also possible [24] and produce qualitatively identical behavior, with a common choice being a simple Poisson sample

$$n_i(x, t + \Delta t) \sim \text{Poisson} \left( \frac{\tilde{n}_i(x, t) \Delta x}{\sigma^2 \Delta t} \right) \frac{\sigma^2 \Delta t}{\Delta x}. \qquad (\text{E11})$$

In the main text, all simulations except for those of Fig. 5 were performed using Eq. (E10). Data for Fig. 5 were generated using Eq. (E11). The final step is to include spatial migrations, which we incorporate using a mass-action term

$$n_i(x, t + \Delta t) \leftarrow n_i(x, t + \Delta t) + \sum_j [K_{ij} n_j(x, t + \Delta t) - K_{ji} n_i(x, t + \Delta t)] \Delta t. \qquad (\text{E12})$$
APPENDIX F: ANALYTICAL TREATMENT
OF DIVERSITY DYNAMICS

The first step to determining the optimal migration rate is to calculate how long on average it takes for an outbreak in deme 1 to spread via migration to deme 2. The model defined by Eqs. (1)–(4) is complex; we will make simplifying assumptions when necessary in order to make progress. We will use the OMD framework of Appendix A to compute the number of infected individuals $N_i$ and antigenic diversity $V_i$ in each deme. We will assume that in deme 1 the viral abundance follows an exponential growth followed by an exponential decay

$$N_1(t) = \begin{cases} N_0 e^{Ft} & \text{if } t \le T_1 \\ N_{\max} e^{-\gamma(t-T_1)} & \text{if } t > T_1. \end{cases} \quad (F1)$$

This says that the infected number in deme 1 will grow exponentially with rate $F = \beta - \gamma$ until reaching a maximum at the outbreak peak time $T_1$. At this time, there are $N_{\max} = N_0 e^{Ft_1}$ infected individuals and the infected population will then decay as individuals recover at rate $\gamma$. In deme 2, migrant infected populations will be subject to growth and noise. Our primary interest is in the emergence of the infected population in deme 2 as it grows to avoid stochastic extinction. We will assume that the small viral populations in deme 2 during the early outbreak elicit negligible immune response, instead enjoying a constant growth rate $F$. Integrating Eq. (1) over $x$ and ignoring migrations, we find

$$\frac{dN_2(t)}{dt} = FN_2 + \sigma \sqrt{N_2} \xi(t), \quad (F2)$$

where $\xi(t)$ is a white noise. Using Eq. (F2), we can compute the probability that a virus with initial abundance $N_2(0)$ eventually goes extinct due to noise. The calculation is fastest using common results from martingale theory [78]. The stochastic process $Z(t) = \exp[-2FN_2(t)/\sigma^2]$ is a martingale which is shown by computing $dZ/dt$ using Itô's lemma and noticing that $Z(t)$ has no deterministic rate of change. The process $Z(t)$ being a martingale implies that its expectation should not change in time (with expectation $\langle \cdot \cdot \cdot \rangle$ here being over noise realizations, not over the antigenic distributions as previously discussed). Because $Z(t)$ has a known deterministic initial condition, we immediately have

$$\langle Z(t) \rangle = Z(0) = P(\text{extinction}) \exp\left(-\frac{2F \times 0}{\sigma^2}\right) + P(\text{establishment}) \exp\left(-\frac{2F \times \infty}{\sigma^2}\right), \quad (F3)$$

which immediately gives the probability of establishment 1 – $P(\text{extinction})$ as

$$P(\text{establishment}) = 1 - Z(0) \approx \frac{2F}{\sigma^2} N_2(0). \quad (F4)$$

To include migration, we note that in a given time interval $t \to t + dt$, $k \times N_1(t) \times dt$ individuals migrate from deme 1 to deme 2. Denoting by $\Delta T$ the time at which the first

successful migration arrives in deme 2, we can explicitly calculate

$$P(\Delta T > t) = \prod_{s<t} \left(1 - \frac{2F}{\sigma^2} k N_1(t) dt\right) \quad (F5)$$

$$= \exp\left(-\int_{0}^{t} \frac{2kF}{\sigma^2} N_1(s) ds\right). \quad (F6)$$

Computing the mean of a random variable distributed according to (F6) is difficult in general. Instead, we compute the most likely value for $\Delta T$, which is the inflection point for the cumulative distribution function $\frac{d^2P(\Delta T > t)}{dt^2} = 0$. The most likely value for the migration time $\Delta T$ is given by

$$N_1(\Delta T)^2 = \frac{\sigma^2}{2kF} N_1'(\Delta T). \quad (F7)$$

Assuming that migration happens before the outbreak peaks in deme 1 ($\Delta T < T_1$), we find a simple expression for the migration time

$$\Delta T = \frac{1}{F} \log \left( \frac{\sigma^2}{2kN_0} \right). \quad (F8)$$

We now have a concrete relationship between the rate of host migration and the time to seed an outbreak in deme 2. One could also compute the expected value of $\Delta T$ numerically, $\langle \Delta T \rangle = \int_0^\infty P(\Delta T > t) dt$. Once the first infected individual is successfully established in deme 2, the growth should be approximately exponential $N_2(t) = e^{F(t-\Delta T)}$.

Taking $N_2$ to be exponential but shifted in time by an amount $\langle \Delta T \rangle$, we can evaluate the integrals in Eq. (B6), finding

$$\begin{aligned} V_2(T_1 + \langle \Delta T \rangle) - V_1(T_1) & \approx 2D \langle \Delta T \rangle \{1 - \exp[-kN_0 e^{\langle \Delta T \rangle} (T_1 - \langle \Delta T \rangle + \beta^{-1})]\} \\ & \approx 2D \langle \Delta T \rangle \left[1 - \exp\left(-\frac{\sigma^2}{2}(T_1 - \langle \Delta T \rangle + \beta^{-1})\right)\right]. \end{aligned} \quad (F9)$$

In the last line we have used Eq. (F8) to simplify. While the above function does not admit a simple expression for the optimal migration time $\Delta T$, the exponential can be expanded as

$$V_2(T_1 + \langle \Delta T \rangle) - V_1(T_1) \approx 2D \Delta T \left( \frac{\sigma^2}{2} (T_1 - \Delta T + \beta^{-1}) \right), \quad (F10)$$

from which we expect the optimal waiting time to be

$$\Delta T^* \approx \frac{T1 + \beta^{-1}}{2}. \quad (F11)$$

Using Eq. (F8), we now have an estimate for the optimal migration rate

$$k^* \approx \frac{2N_0}{\sigma^2} \exp\left(F \frac{T1 + \beta^{-1}}{2}\right). \quad (F12)$$

For the parameters used in Fig. 3, this expression gives $k^* \sim 4 \times 10^{-5}$.
## APPENDIX G: MIGRATION RATE PREDICTIONS

In Fig. 7(b) we compare migration rates between different regions in the world to see what probability of viral survival we would have in these systems. To do this, we obtain data on the number of flights per year within that system (for the U.S. calculation we would find the number of internal U.S. flights per year) as well as the total number of people in that system (number of people living in the U.S.) and the average number of passengers on a plane flying within the region. The calculation would follow as

$$ migration\ rate = \frac{(\text{No. of flights per year})(\text{average No. of people per flight})}{(52\text{ weeks per year})(\text{No. of people in region})} $$

We end up with the number of people moving per week divided by the total people in the region as a migration rate prediction. The number of people traveling per week is used due to the migration rate scaling with the recovery time, which is about a week. It is important to acknowledge the approximative nature of these calculations as the people traveling may not necessarily live in the place the flight originates from and therefore not belong to the number of people in that region.

## APPENDIX H: U.S. MIGRATION RATE CALCULATION

From the following calculation, travel within the U.S. is well mixed:

$$ U.S.\ migration\ rate = \frac{(16\ 405\ 000\ flights\ per\ year)(200\ passenger\ per\ flight)}{(52\ weeks\ per\ year)(333\ 300\ 000\ U.S.\ residents)} = 2 \times 10^{-1}. $$

## APPENDIX I: MIGRATION RATE AND POPULATION SIZE HETEROGENEITY

We expand the airport network analysis to the top 30 most trafficked airports. To determine whether the specific network topology is necessary for the observation of the negative trend, we conduct two sets of simulations following the same procedure: We (i) combine all migration rates into one distribution, (ii) randomly sample with replacement this migration distribution for all entries in the adjacency matrix, (iii) ensure all diagonal elements are zero, (iv) set all initial population sizes for each of the 30 demes to the same $N_h$ value, and (v) find the survival probability for the virus for an outbreak which occurs at each deme in the network.

When the assumption that the flux into an airport is roughly the same as the flux out is enforced, the adjacency matrix is always symmetric. To test if this structure is imperative to produce the negative relationship, we randomly sample half of the migration rates from the distribution and reflect these rates across the diagonal to establish symmetry. In Fig. 8(d) the relationship between the closeness centrality and the survival probability is restored, indicating that what is important is that the fluxes in and out of each deme are roughly equal.

[1] G. L. Kaján, A. Doszpoly, Z. L. Tarján, M. Z. Vidovszky, and T. Papp, Virus-host coevolution with a focus on animal and human DNA viruses, J. Mol. Evol. 88, 41 (2020).

[2] L. M. Sompayrac, How the Immune System Works (Wiley, New York, 2022).

[3] R. Ahmed and D. Gray, Immunological memory and protective immunity: Understanding their relation, Science 272, 54 (1996).

[4] C. O. McCoy, T. Bedford, V. N. Minin, P. Bradley, H. Robins, and F. A. Matsen IV, Quantifying evolutionary constraints on B-cell affinity maturation, Philos. Trans. Roy. Soc. B 370, 20140244 (2015).

[5] A. M. Phillips, K. R. Lawrence, A. Moulana, T. Dupic, J. Chang, M. S. Johnson, I. Cvijovic, T. Mora, A. M. Walczak, and M. M. Desai, Binding affinity landscapes constrain the evolution of broadly neutralizing anti-influenza antibodies, eLife 10, e71393 (2021).

[6] A. Mazzolini, T. Mora, and A. M. Walczak, Inspecting the interaction between human immunodeficiency virus and the immune system through genetic turnover, Philos. Trans. R. Soc. B 378, 20220056 (2023).

[7] V. Chardès, M. Vergassola, A. M. Walczak, and T. Mora, Affinity maturation for an optimal balance between long-

term immune coverage and short-term resource constraints, Proc. Natl. Acad. Sci. USA 119, e2113512119 (2022).

[8] B. Thyagarajan and J. D. Bloom, The inherent mutational tolerance and antigenic evolvability of influenza hemagglutinin, eLife 3, e03300 (2014).

[9] D. H. Morris, K. M. Gostic, S. Pompei, T. Bedford, M. Łuksza, R. A. Neher, B. T. Grenfell, M. Lässig, and J. W. McCauley, Predictive modeling of influenza shows the promise of applied evolutionary biology, Trends Microbiol. 26, 102 (2018).

[10] T. Sakai and Y. Morimoto, The history of infectious diseases and medicine, Pathogens 11, 1147 (2022).

[11] M. Łuksza and M. Lässig, A predictive fitness model for influenza, Nature (London) 507, 57 (2014).

[12] T. Bedford, S. Cobey, P. Beerli, and M. Pascual, Global migration dynamics underlie evolution and persistence of human influenza a (H3N2), PLoS Pathog. 6, e1000918 (2010).

[13] N. Strelkowa and M. Lässig, Clonal interference in the evolution of influenza, Genetics 192, 671 (2012).

[14] J. D. Bloom, L. I. Gong, and D. Baltimore, Permissive secondary mutations enable the evolution of influenza oseltamivir resistance, Science 328, 1272 (2010).

[15] M. B. Doud, J. M. Lee, and J. D. Bloom, How single mutations affect viral escape from broad and narrow
antibodies to H1 influenza hemagglutinin, Nat. Commun. 9, 1386 (2018).

[16] J. M. Lee, J. Huddleston, M. B. Doud, K. A. Hooper, N. C. Wu, T. Bedford, and J. D. Bloom, Deep mutational scanning of hemagglutinin helps predict evolutionary fates of human H3N2 influenza variants, Proc. Natl. Acad. Sci. USA 115, E8276 (2018).

[17] V. N. Petrova and C. A. Russell, The evolution of seasonal influenza viruses, Nat. Rev. Microbiol. 16, 47 (2018).

[18] D. J. Smith, Mapping the antigenic and genetic evolution of influenza virus, Science 305, 371 (2004).

[19] T. Bedford, A. Rambaut, and M. Pascual, Canalization of the evolutionary trajectory of the human influenza virus, BMC Biol. 10, 38 (2012).

[20] A. Sasaki, Evolution of antigen drift/switching: Continuously evading pathogens, J. Theor. Biol. 168, 291 (1994).

[21] L. Yan, R. Neher, and B. I. Shraiman, Phylodynamics of rapidly adapting pathogens: Extinction and speciation of a Red Queen, arXiv:1810.11918.

[22] J. Marchi, M. Lässig, T. Mora, and A. M. Walczak, Multi-lineage evolution in viral populations driven by host immune systems, Pathogens 8, 115 (2019).

[23] J. Marchi, M. Lässig, A. M. Walczak, and T. Mora, Antigenic waves of virus-immune coevolution, Proc. Natl. Acad. Sci. USA 118, e2103398118 (2021).

[24] V. Chardès, A. Mazzolini, T. Mora, and A. M. Walczak, Evolutionary stability of antigenically escaping viruses, Proc. Natl. Acad. Sci. USA 120, e2307712120 (2023).

[25] D. Balcan, V. Colizza, B. Gonçalves, H. Hu, J. J. Ramasco, and A. Vespignani, Multiscale mobility networks and the spatial spreading of infectious diseases, Proc. Natl. Acad. Sci. USA 106, 21484 (2009).

[26] G. De Luca, K. Van Kerchove, P. Coletti, C. Poletto, N. Bossuyt, N. Hens, and V. Colizza, The impact of regular school closure on seasonal influenza epidemics: A data-driven spatial transmission model for Belgium, BMC Infect. Dis. 18, 29 (2018).

[27] G. Pullano et al., Underdetection of cases of COVID-19 in France threatens epidemic control, Nature (London) 590, 134 (2021).

[28] C. Poletto, S. Meloni, V. Colizza, Y. Moreno, and A. Vespignani, Host mobility drives pathogen competition in spatially structured populations, PLoS Comput. Biol. 9, e1003169 (2013).

[29] T. Bedford et al., Global circulation patterns of seasonal influenza viruses vary with antigenic drift, Nature (London) 523, 217 (2015).

[30] F. Wen, T. Bedford, and S. Cobey, Explaining the geographical origins of seasonal influenza A (H3N2), Proc. R. Soc. B 283, 20161312 (2016).

[31] J. Hadfield, C. Megill, S. M. Bell, J. Huddleston, B. Potter, C. Callender, P. Sagulenko, T. Bedford, and R. A. Neher, Nextstrain: Real-time tracking of pathogen evolution, Bioinformatics 34, 4121 (2018).

[32] J. L.-H. Tsui et al., Genomic assessment of invasion dynamics of SARS-CoV-2 Omicron BA.1, Science 381, 336 (2023).

[33] M. J. Keeling, Metapopulation moments: Coupling, stochasticity and persistence, J. Anim. Ecol. 69, 725 (2000).

[34] J. W. Fox, D. Vasseur, M. Cotroneo, L. Guan, and F. Simon, Population extinctions can increase metapopulation persistence, Nat. Ecol. Evol. 1, 1271 (2017).

[35] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Epidemic processes in complex networks, Rev. Mod. Phys. 87, 925 (2015).

[36] M. F. Boni, J. R. Gog, V. Andreasen, and M. W. Feldman, Epidemic dynamics and antigenic evolution in a single season of influenza A, Proc. R. Soc. B 273, 1307 (2006).

[37] R. Kumata and A. Sasaki, Antigenic escape is accelerated by the presence of immunocompromised hosts, Proc. R. Soc. B 289, 20221437 (2022).

[38] A. Sasaki, S. Lion, and M. Boots, Antigenic escape selects for the evolution of higher pathogen transmission and virulence, Nat. Ecol. Evol. 6, 51 (2022).

[39] M. Jesse and H. Heesterbeek, Divide and conquer? Persistence of infectious agents in spatial metapopulations of hosts, J. Theor. Biol. 275, 12 (2011).

[40] L. C. McManus et al., Evolution reverses the effect of network structure on metapopulation persistence, Ecology 102, e03381 (2021).

[41] A. Lamarins, S. M. Carlson, E. Prévost, W. H. Satterthwaite, and M. Buoro, Eco-evolutionary consequences of selective exploitation on metapopulations illustrated with Atlantic salmon, Fish Fish. 26, 291 (2025).

[42] B. Houchmandzadeh and M. Vallade, The fixation probability of a beneficial mutation in a geographically structured population, New J. Phys. 13, 073020 (2011).

[43] M. Slatkin, Populational heritability, Evolution 35, 859 (1981).

[44] L. Marrec, I. Lamberti, and A. F. Bitbol, Toward a universal model for spatially structured populations, Phys. Rev. Lett. 127, 218102 (2021).

[45] Y. P. Kuo and O. Carja, Evolutionary graph theory beyond single mutation dynamics: On how network-structured populations cross fitness landscapes, Genetics 227, iyae055 (2024).

[46] L. A. Segel and A. S. Perelson, Shape space: An approach to the evaluation of cross-reactivity effects, stability and controllability in the immune system, Immunol. Lett. 22, 91 (1989).

[47] R. M. May and R. M. Anderson, Population biology of infectious diseases: Part II, Nature (London) 280, 455 (1979).

[48] R. M. Anderson and R. M. May, Population biology of infectious diseases: Part I, Nature (London) 280, 361 (1979).

[49] R. M. Anderson and R. M. May, Directly transmitted infections diseases: Control by vaccination, Science 215, 1053 (1982).

[50] R. Anderson, *Infectious Diseases of Humans: Dynamics and Control* (Cambridge University Press, Cambridge, 1991).

[51] WorldPop, Global Flight Data Annual, https://hub.worldpop.org/project/categories?id=13 (2015).

[52] Statistics and market data, https://www.statista.com (Statista, New York, 2024).

[53] Global metrics: Cities, https://www.macrotrends.net/global-metrics/cities (Macrotrends, Seattle, 2024).

[54] M. T. Pearce, A. Agarwala, and D. S. Fisher, Stabilization of extensive fine-scale diversity by ecologically driven spatiotemporal chaos, Proc. Natl. Acad. Sci. USA 117, 14572 (2020).

[55] S. Gandon, Y. Capowiez, Y. Dubois, Y. Michalakis, and I. Olivieri, Local adaptation and gene-for-gene coevolution in a metapopulation model, Proc. R. Soc. B 263, 1003 (1996).

[56] A. Aleta, A. N. S. Hisi, S. Meloni, C. Poletto, V. Colizza, and Y. Moreno, Human mobility networks and persistence of rapidly mutating pathogens, R. Soc. Open Sci. 4, 160914 (2017).
[57] A. N. S. Hisi, E. E. N. Macau, and L. H. G. Tizei, The role of mobility in epidemic dynamics, Physica A 526, 120663 (2019).

[58] P. L. Thompson and E. A. Fronhofer, The conflict between adaptation and dispersal for maintaining biodiversity in changing environments, Proc. Natl. Acad. Sci. USA 116, 21061 (2019).

[59] S. R. Meakin and M. J. Keeling, Correlations between stochastic epidemics in two interacting populations, Epidemics 26, 58 (2019).

[60] M. J. Keeling and P. Rohani, Estimating spatial coupling in epidemiological systems: A mechanistic approach, Ecol. Lett. 5, 20 (2002).

[61] M. Pautasso, M. Moslonka-Lefebvre, and M. J. Jeger, The number of links to and from the starting node as a predictor of epidemic size in small-size directed networks, Ecol. Complex. 7, 424 (2010).

[62] D. H. Gent, S. Bhattacharyya, and T. Ruiz, Prediction of spread and regional development of hop powdery mildew: A network analysis, Phytopathology 109, 1392 (2019).

[63] A. N. Bishop and I. Shames, Link operations for slowing the spread of disease in complex networks, Europhys. Lett. 95, 18005 (2011).

[64] A. Darabi and M. Siami, Centrality-based traffic restriction in delayed epidemic networks, SIAM J. Appl. Dyn. Syst. 22, 3165 (2023).

[65] S. Chang, E. Pierson, P. W. Koh, J. Gerardin, B. Redbird, D. Grusky, and J. Leskovec, Mobility network models of COVID-19 explain inequities and inform reopening, Nature (London) 589, 82 (2021).

[66] K. G. C. Smith, A. Light, G. J. V. Nossal, and D. M. Tarlinton, The extent of affinity maturation differs between the memory and antibody-forming cell compartments in the primary immune response, EMBO J. 16, 2996 (1997).

[67] D. J. Smith, S. Forrest, R. R. Hightower, and A. S. Perelson, Deriving shape space parameters from immunological data, J. Theor. Biol. 189, 141 (1997).

[68] T. Bedford, M. A. Suchard, P. Lemey, G. Dudas, V. Gregory, A. J. Hay, J. W. McCauley, C. A. Russell, D. J. Smith, and A. Rambaut, Integrating influenza antigenic dynamics with molecular evolution, eLife 3, e01914 (2014).

[69] J. Moore, H. Ahmed, and R. Antia, High dimensional random walks can appear low dimensional: Application to influenza H3N2 evolution, J. Theor. Biol. 447, 56 (2018).

[70] https://github.com/caelan-brooks/viral-coev-net

[71] S. Lion, A. Sasaki, and M. Boots, Extending eco-evolutionary theory with oligomorphic dynamics, Ecol. Lett. 26, S22 (2023).

[72] A. Sasaki and U. Dieckmann, Oligomorphic dynamics for analyzing the quantitative genetics of adaptive speciation, J. Math. Biol. 63, 601 (2011).

[73] L. Pechenik and H. Levine, Interfacial velocity corrections due to multiplicative noise, Phys. Rev. E 59, 3893 (1999).

[74] H. Weissmann, N. M. Shnerb, and D. A. Kessler, Simulation of spatial systems with demographic noise, Phys. Rev. E 98, 022131 (2018).

[75] I. Dornic, H. Chaté, and M. A. Munoz, Integration of Langevin equations with multiplicative noise and the viability of field theories for absorbing phase transitions, Phys. Rev. Lett. 94, 100601 (2005).

[76] R. I. Issa, Solution of the implicitly discretised fluid flow equations by operator-splitting, J. Comput. Phys. 62, 40 (1986).

[77] S. MacNamara and G. Strang, in Splitting Methods in Communication, Imaging, Science, and Engineering, edited by R. Glowinski, S. J. Osher, and W. Yin (Springer, Cham, 2016), pp. 95–114.

[78] É. Roldán, I. Neri, R. Chetrite, S. Gupta, S. Pigolotti, F. Jülicher, and K. Sekimoto, Martingales for physicists: A treatise on stochastic thermodynamics and beyond, Adv. Phys. 72, 1 (2023).