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

**HAL**
open science

# The speed of vaccination rollout and the risk of pathogen adaptation

Sylvain Gandon, Amaury Lambert, Troy Day, Todd L Parsons

► To cite this version:

Sylvain Gandon, Amaury Lambert, Troy Day, Todd L Parsons. The speed of vaccination rollout and the risk of pathogen adaptation. 2022. hal-03849469

HAL Id: hal-03849469

https://cnrs.hal.science/hal-03849469

Preprint submitted on 11 Nov 2022

**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.
# The speed of vaccination rollout and the risk of pathogen adaptation

Sylvain Gandon, Amaury Lambert, Marina Voinson, Troy Day, Todd L. Parsons

► To cite this version:

Sylvain Gandon, Amaury Lambert, Marina Voinson, Troy Day, Todd L. Parsons. The speed of vaccination rollout and the risk of pathogen adaptation. Journal of the Royal Society Interface, 2025, 22 (228), pp.20250060. 10.1101/2022.08.01.22278283. hal-03849469v2

HAL Id: hal-03849469

https://cnrs.hal.science/hal-03849469v2

Submitted on 3 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
# The speed of vaccination rollout and the risk of pathogen adaptation

Sylvain Gandon¹*, Amaury Lambert²,³, Troy Day⁴,⁵, and Todd L. Parsons⁶

¹CEFE, CNRS, Univ Montpellier, EPHE, IRD, Montpellier, France

*Corresponding author

²Institut de Biologie de l'ENS (IBENS), École Normale Supérieure (ENS), CNRS UMR 8197, Paris, France

³Center for Interdisciplinary Research in Biology (CIRB), Collège de France, CNRS UMR 7241, INSERM U1050, PSL Research University, Paris, France

⁴Department of Mathematics and Statistics, Queen's University, Kingston, Canada

⁵Department of Biology, Queen's University, Kingston, Canada

⁶Laboratoire de Probabilités, Statistique et Modélisation (LPSM), Sorbonne Université, CNRS UMR 8001, Paris, France

Saturday 13th August, 2022

E-mails: SG - sylvain.gandon@cefe.cnrs.fr; AL - amaury.lambert@ens.fr; TD - tday@mast.queensu.ca; TLP - todd.parsons@upmc.fr;

Keywords: evolutionary epidemiology, vaccination, life-history evolution, demographic stochasticity, virulence, adaptive dynamics
INTERFACE

royalsocietypublishing.org/journal/rsif

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

Research

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

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

Check for
updates

**Cite this article:** Gandon S, Lambert A, Voinson M, Day T, Parsons TL. 2025 The speed of vaccination rollout and the risk of pathogen adaptation. *J. R. Soc. Interface* **22**: 20250060.

https://doi.org/10.1098/rsif.2025.0060

Received: 17 January 2025

Accepted: 28 April 2025

**Subject Category:**
Life Sciences–Mathematics interface

**Subject Areas:**
evolution

**Keywords:**
vaccination, evolutionary epidemiology, demographic stochasticity, adaptive dynamics

**Author for correspondence:**
Sylvain Gandon
e-mail: sylvain.gandon@cefe.cnrs.fr

Electronic supplementary material is available
online at https://doi.org/10.6084/
m9 figshare.c.7828935.

# The speed of vaccination rollout and the risk of pathogen adaptation

Sylvain Gandon<fnref n="1" />, Amaury Lambert<fnref n="2" />,<fnref n="3" />, Marina Voinson<fnref n="1" />, Troy Day<fnref n="4" />,⁵ and Todd L. Parsons<fnref n="6" />

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

SG, 0000-0003-2624-7856; AL, 0000-0002-7248-9955; MV, 0000-0001-6085-9974; TD, 0000-0002-1052-6140; TLP, 0000-0002-2599-8415

Vaccination is expected to reduce disease prevalence and to halt the spread of epidemics. Pathogen adaptation, however, may erode the efficacy of vaccination and limit our ability to control disease spread. Here, we examine the influence of the rate of vaccination of the host population on the overall risk of pathogen adaptation to vaccination. We extend the framework of evolutionary epidemiology theory to account for demographic stochasticity in the different steps leading to the adaptation to vaccination: (i) the introduction of a vaccine-escape variant by mutation from an endemic wild-type pathogen, (ii) the invasion of this vaccine-escape variant in spite of the risk of early extinction, (iii) the spread and fixation of the vaccine-escape variant in the pathogen population. We introduce a novel and versatile hybrid analytical-numerical method that allows fast computation of the probabilities associated with these steps. Using it, we show that increasing the rate of vaccination can reduce both the number of cases and the likelihood of pathogen adaptation. Our work clarifies the influence of vaccination policies—a major ecological perturbation of the environment of a pathogen—on different steps of pathogen adaptation. The model provides a useful theoretical framework to account for the interplay between epidemiology, selection and genetic drift and to anticipate the effects of public-health interventions on pathogen evolution.

## 1. Introduction

Vaccination offers unique opportunities to protect a large fraction of the host population and thus to control spreading epidemics. In principle, comprehensive vaccination coverage can lead to pathogen eradication. In practice, however, the coverage required for eradication is often impossible to reach with imperfect vaccines [1,2]. Moreover, pathogen adaptation may erode the efficacy of vaccination. Even if adaptation to vaccines is less common than adaptation to drugs [3–5] the spread of vaccine-escape mutations may challenge our ability to halt the spread of epidemics.

Understanding the dynamics of pathogen adaptation to vaccines is particularly relevant in the control of many human and animal pathogens [6–10]. Yet, most theoretical studies that explore the evolution of pathogens after vaccination are based on the analysis of deterministic models and ignore the potential effects induced by the stochasticity of epidemiological dynamics [7,11,12].

© 2025 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

<footnote n="1">¹CEFE, CNRS, Univ Montpellier, EPHE, IRD, Montpellier, France</footnote>

<footnote n="2">²Institut de Biologie de l'ENS (IBENS), Ecole Normale Supérieure (ENS), CNRS UMR 8197, PSL Research University, Paris, France</footnote>

<footnote n="3">³Center for Interdisciplinary Research in Biology (CIRB), Collège de France, CNRS UMR 7241, INSERM U1050, PSL Research University, Paris, France</footnote>

<footnote n="4">⁴Department of Mathematics and Statistics, and ⁵Department of Biology, Queen's University, Kingston, Ontario, Canada</footnote>

<footnote n="6">⁶Sorbonne Université, Université Paris Cité, CNRS, Laboratoire de Probabilités, Statistique et Modélisation, LPSM, Paris, France</footnote>
Abstract

Vaccination is expected to reduce disease prevalence and to halt the spread of epidemics. But pathogen adaptation may erode the efficacy of vaccination and challenge our ability to control disease spread. Here we examine the influence of the speed of vaccination rollout on the overall risk of pathogen adaptation to vaccination. We extend the framework of evolutionary epidemiology theory to account for the different steps leading to adaptation to vaccines: (1) introduction of a vaccine-escape variant by mutation from an endemic wild-type pathogen, (2) invasion of this vaccine-escape variant in spite of the risk of early extinction, (3) spread and, eventually, fixation of the vaccine-escape variant in the pathogen population. We show that the risk of pathogen adaptation is maximal for intermediate speed of vaccination rollout. On the one hand, slower rollout decreases pathogen adaptation because selection is too weak to avoid early extinction of the new variant. On the other hand, faster rollout decreases pathogen adaptation because it reduces the influx of adaptive mutations. Hence, vaccinating faster is recommended to decrease both the number of cases and the likelihood of pathogen adaptation. We also show that pathogen adaptation is driven by its basic reproduction ratio, the efficacy of the vaccine and the effects of the vaccine-escape mutations on pathogen life-history traits. Accounting for the interplay between epidemiology, selection and genetic drift, our work clarifies the influence of vaccination policies on different steps of pathogen adaptation and allows us to anticipate the effects of public-health interventions on pathogen evolution.

**Significance statement:** Pathogen adaptation to host immunity challenges the efficacy of vaccination against infectious diseases. Are there vaccination strategies that limit the emergence and the spread of vaccine-escape variants? Our theoretical model clarifies the interplay between the timing of vaccine escape mutation events and the transient epidemiological dynamics following the start of a vaccination campaign on pathogen adaptation. We show that the risk of adaptation is maximized for intermediate vaccination coverage but can be reduced by a combination of non pharmaceutical interventions and maximizing the speed of the vaccination rollout. These recommendations may have important implications for the choice of vaccination strategies against the ongoing SARS-CoV-2 pandemic.
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Graphical representation of the different evolutionary epidemiology outcomes after vaccination. The density of the wild-type pathogen is indicated in light blue and the dynamics of the mutant in orange. Each panel describes the temporal dynamics of the epidemics after the start of vaccination: (a) eradication of the wild-type pathogen, (b) new endemic equilibrium of the wild-type population after damped oscillations (with no introduction of the vaccine-escape mutant), (c) early extinction of the vaccine-escape mutant after its introduction by mutation, (d) invasion of the vaccine-escape mutant followed by its extinction, (e) invasion of the vaccine-escape mutant and long-term coexistence with the wild-type in a new endemic equilibrium after damped oscillations, (f) invasion and fixation of the vaccine-escape mutant (extinction of the wild-type). The vertical dashed line (black) indicates the start of vaccination. For simplicity, we consider that vaccination starts after the wild-type population has reached an endemic equilibrium. The horizontal dashed line indicates the 'stochastic threshold' above which one may consider that the deterministic model provides a very good approximation of the dynamics and we can neglect the effect of demographic stochasticity. *Invasion* occurs when the vaccine-escape variant manages to go beyond the 'stochastic threshold' (d–f). *Adaptation* occurs when the vaccine-escape variant is maintained in the population (e,f). *Fixation* occurs when the vaccine-escape variant manages to outcompete the wild-type (f).</figcaption></figure>

Demographic stochasticity, however, can affect the establishment of new mutations and the long-term evolution of pathogens [6,13,14]. Several studies have shown how the demographic stochasticity induced by finite host and pathogen population sizes alters selection on the life-history traits of pathogens [15–17]. These analytical predictions rely on the assumption that the rate of pathogen mutation is low, which allows one to decouple epidemiological and evolutionary time scales. Indeed, when the influx of new mutations is low, the new strain is always introduced after the resident pathogen population has reached an endemic equilibrium. Many pathogens, however, have relatively high mutation rates [18] and the fate of a pathogen mutant introduced away from the endemic equilibrium is likely to be affected by the dynamics of the pathogen populations. Moreover, the start of a vaccination campaign is expected to yield large perturbations of the epidemiological dynamics and new mutations are likely to appear when the pathogen population is far from its endemic equilibrium.

The aim of the present study is to develop a versatile theoretical framework to evaluate the consequences of vaccination on the risk of pathogen adaptation to vaccination. Vaccination, because it reduces the availability of fully susceptible hosts, constitutes a major perturbation of the environment of the pathogen. There are six main evolutionary-epidemiological outcomes after the start of vaccination which are summarized in figure 1. Some of these outcomes are more desirable because they do not lead to the invasion of a new variant (figure 1a–c). The aim of the present study is to help anticipate what factors may promote the inva- sion of vaccine-escape variants (figure 1e,f). In the following, we use a combination of deterministic and birth-and-death process approximations to study the joint epidemiological and evolutionary dynamics of the pathogen population. This analysis extends the evolutionary epidemiology theory framework [19,20] to account for the stochasticity of the mutation process and the impact of demographic stochasticity on the risk of extinction of a new variant. We explore a relatively simple vaccination model to reveal the importance of the speed of the vaccination rollout as well as of the life-history characteristics of the vaccine-escape variants on three major steps in the adaptation of the pathogen: (i) introduction of a variant by mutation, (ii) invasion and (iii) fixation. More complex epidemiological scenarios are explored in the electronic supplementary material and in §4.

## 2. The vaccination model

We use a classical SIR epidemiological model with vital dynamics (i.e. host births and deaths) [21], where hosts can be susceptible, infected or recovered [22], and are either vaccinated or unvaccinated. A host may be infected by one of two strains: a resident wild-type, or a novel mutant (we assume co-infections are not possible).

We consider a continuous-time Markov process tracking the *number* of individuals of each type of host (see table 1 for a detailed description). Rates are interpreted as probabilities per unit time. We incorporate vital dynamics by assuming that all hosts have a base mortality rate of $\delta$, while new susceptible hosts are recruited at rate $\nu n$. Here, $n$ is a ‘system size’, or scaling parameter, that
21

# 1 Introduction

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

Vaccination offers unique opportunities to protect a large fraction of the host population and thus to control spreading epidemics. In principle, large vaccination coverage can lead to pathogen eradication. In practice, however, the coverage required for eradication is often impossible to reach with imperfect vaccines [16, 34]. Moreover, pathogen adaptation may erode the efficacy of vaccination. Even if adaptation to vaccines is less common than adaptation to drugs [14, 25, 26] the spread of vaccine-escape mutations may challenge our ability to halt the spread of epidemics.

Understanding the dynamics of pathogen adaptation to vaccines is particularly relevant in the control of the ongoing SARS-CoV-2 pandemic. Yet, most theoretical studies that explore the evolution of pathogens after vaccination are based on the analysis of deterministic models and ignore the potential effects induced by the stochasticity of epidemiological dynamics. Demographic stochasticity, however, drives the intensity of genetic drift and can affect the establishment of new mutations and the long-term evolution of pathogens [42, 44, 41]. Several studies showed how the demographic stochasticity induced by finite host and pathogen population sizes alters selection on the life-history traits of pathogens [28, 22, 36]. These analytical predictions rely on the assumption that mutation rate is low, which allows us to decouple epidemiological and evolutionary time scales. Indeed, when mutation rate is low, the new strain is always introduced after the resident pathogen population has reached its endemic equilibrium. Many pathogens, however, have relatively large mutation rates [43] and the fate of a pathogen mutant introduced away from the endemic equilibrium is likely to be affected by the dynamics of the pathogen populations. Besides, the start of a vaccination campaign is expected to yield massive perturbations of the epidemiological dynamics and new mutations are likely to appear when the pathogen population is far from its endemic equilibrium.

The aim of the present study is to develop a versatile theoretical framework to evaluate the consequences of vaccination on the risk of pathogen adaptation to vaccination. There are six main evolutionary epidemiology outcomes after the start of vaccination which are summarized in **Figure 1**. Some of these outcomes are more favorable than others because they do not lead to the invasion of a new variant (**Figure 1a-c**). In contrast, vaccination may lead to the invasion of
<figure><img src="image_7.png" /><figcaption>**Figure 2.** A schematic representation of the model. Naive and uninfected hosts ($S_u^n$ and $S_v^n$) are introduced at a rate $\nu$ and are vaccinated at rate $\nu$. Immunization induced by the vaccine wanes at rate $\omega_v$. Uninfected hosts ($S_u^n$ and $S_v^n$) die at a rate $\delta$ while infected hosts ($I_{ui}^n$ and $I_{vi}^n$) die at a rate $d_i = \delta + \alpha_i$, where $i$ refers to the virus genotype: the wild-type ($i = w$) or the vaccine-escape mutant ($i = m$). The rate of infection of naive hosts by the genotype $i$ is $h_i = \beta_i(I_{ui}^n + I_{vi}^n)/H^n$, where $\beta_i$ is the transmission rate of the genotype $i$. Vaccination reduces the force of infection and $\epsilon_i$ refers to the ability of the genotype $i$ to escape the immunity triggered by vaccination (we assume $\epsilon_m > \epsilon_w$). A host infected by pathogen genotype $i$ recovers from the infection at rate $\gamma_i$ and yields naturally immune hosts ($R^n$ hosts) that cannot be reinfected by both the wild-type and the escape mutant. Natural immunity is assumed to wane at rate $\omega_r$. The total host population density is $H^n = S_u^n + S_v^n + \sum_{i \in \{w, m\}} (I_{ui}^n + I_{vi}^n) + R^n$.</figcaption></figure>

indicates the order of magnitude of the arena in which the epidemic occurs: the total host population size varies stochastically in time, but remains of the order of $n$ individuals. We track the numbers of two classes of susceptible hosts, unvaccinated, u, or vaccinated, v, ($S_u^n$, $S_v^n$), four classes of unvaccinated and vaccinated individuals, infected with the wild-type, w, ($I_{uw}^n$, $I_{vw}^n$) or with a mutant strain, m, ($I_{um}^n$, $I_{vm}^n$), and the number recovered, $R^n$. The total number susceptible is thus $S^n = S_u^n + S_v^n$, while the number of infected hosts is $I^n = \sum_{i \in \{w, m\}} I_{ui}^n + I_{vi}^n$. We write $H^n$ for the total number of hosts:

$$H^n = S^n + I^n + R^n. \qquad (2.1)$$

Vaccination is assumed to take place at a constant rate $\nu$ for all susceptible hosts. The immunity triggered by vaccination is assumed to wane at rate $\omega_v$, and natural (i.e. infection-induced) immunity is assumed to wane at rate $\omega_r$. We further assume that the virulence $\alpha_i$ (the mortality rate induced by the infection), the transmission $\beta_i$ (the production rate of new infections) and the recovery $\gamma_i$ (the rate at which the host clears the infection) are fully governed by the pathogen genotype ($i = w$ or $m$). A fourth trait, $\epsilon_i \in [0, 1]$, governs the infectivity of pathogen genotype $i$ on vaccinated hosts (infectivity of all genotypes is assumed to be equal to 1 on unvaccinated hosts). In other words, this final trait measures the ability of the pathogen to escape the immunity triggered by the vaccine, but once a pathogen infects a vaccinated hosts, its phenotype (transmission, virulence and recovery rates) is assumed to be the same as in naive infected hosts. This assumption allows us to aggregate infected hosts irrespective of their vaccination status, which simplifies the analysis below: when there is a single class of infected hosts it is easier to describe how epidemiological dynamics feed back onto pathogen adaptation. Our model, however, could be readily extended to more realistic scenarios where, for instance, transmission and virulence are reduced in vaccinated hosts (see §4). In the simplest version of our model, we assume that recovered individuals are fully protected (no reinfections are allowed before immunity wanes) because natural immunity is expected to be more effective than immunity triggered by vaccination (e.g. this is believed to be true for measles [24] and influenza [25–27], but not necessarily for SARS-CoV-2 [28]). Again, we make this assumption for the sake of simplicity, but we relax it in the electronic supplementary material, where we also discuss the robustness of our results when vaccine-escape mutants are allowed to reinfect naturally immune hosts. We further assume frequency-dependent transmission where the number of contacts a host may have in the population is constant, but a proportion of those contacts may be infectious. Note, however, that other forms of transmission (e.g. density-dependent transmission where the number of contacts a host may have in the population is constant, but a proportion of those contacts may be infectious [29]) are expected to affect quantitatively the spread of the epidemic but to yield qualitatively similar results on pathogen evolution [30]. We summarize the states of the process and the jump rates at which individuals transition between states in table 1 and in figure 2.

We use this model to examine the epidemiological and evolutionary dynamics following the start of a vaccination campaign. For the sake of simplicity, we focus our analysis on scenarios where the pathogen population has reached an endemic equilibrium before the start of vaccination. This is a strong assumption, but our aim in this study is to focus on a scenario where the initial epidemiological state of the system is fixed to understand the stochastic fate of vaccine-escape mutations during the transient epidemiological dynamics of the pathogen population following the start of the vaccination campaign. This is a necessary first step before studying more complex scenarios where vaccination starts before the epidemic has reached an endemic equilibrium (a scenario to which our approach is readily adapted). The default parameter values used to explore numerically the dynamics of
vaccine-escape variants (Figure 1e-f). In the following we use a combination of deterministic and branching process approximations to study the joint epidemiological and evolutionary dynamics of the pathogen population. This analysis reveals the importance of the speed of the vaccination rollout as well as of the life-history characteristics of the vaccine-escape variants on the probability of pathogen adaptation.

## 2 Model

We use a classical SIR epidemiological model with demography, where hosts can either be susceptible, infected or recovered [3]. The discrete number of each of these types of hosts is denoted by $S^n$, $I^n$ and $R^n$, respectively. Because we are interested in the effect of demographic stochasticity the model is derived from a microscopic description of all the events that may occur in a finite—but not fixed—host population of (varying) total size $N^n = S^n + I^n + R^n$. We consider a continuous-time Markov process tracking the number of individuals of each type of host (see SI Section 1 for a detailed description). The susceptible hosts immigrate at rate $\lambda n$, where n is a “system size”, or scaling parameter, that indicates the order of magnitude of the arena in which the epidemic occurs. Hence the total host population varies stochastically in time, but remains of the order of n.

Vaccination may either take place with probability p when a new susceptible host enters the population (e.g., vaccination of newborns) or at a constant rate v for all other susceptible hosts (e.g., vaccination of adults). The immunity triggered by vaccination is assumed to wane at rate ω. A host may be unvaccinated, U, or vaccinated, V, and may either be uninfected or infected with the wild type, w, or a mutant strain, m (we assume coinfections are not possible). We thus have to track the numbers of two classes of susceptible hosts ($S_U^n, S_V^n$) and four classes of infected individuals ($I_{Uw}^n, I_{Um}^n, I_{Vw}^n, I_{Vm}^n$). We assume that the virulence $\alpha_i$ (the mortality rate induced by the infection), the transmission $\beta_i$ (the production rate of new infections), and the recovery $\gamma_i$ (the rate at which the host clears the infection) are fully governed by the pathogen genotype (i = w or m). A fourth trait $e_i \in [0, 1]$ governs the infectivity of pathogen genotype i on vaccinated hosts (infectivity of all genotypes is assumed to be equal to 1 on unvaccinated hosts). In other words, this final trait measures the ability of the pathogen to escape the immunity triggered by the vaccine.
<table><caption>**Table 1.** We model the epidemic via a continuous time Markov chain (CTMC) with discrete states ($S_u^n, S_v^n, I_{uw}^n, I_{um}^n, I_{vw}^n, I_{vm}^n, R^n$). Jumps ($\Delta S_u^n, \Delta S_v^n, \Delta I_{uw}^n, \Delta I_{um}^n, \Delta I_{vw}^n, \Delta I_{vm}^n, \Delta R^n$) occur at state dependent rates (i.e. with probability proportional to $\Delta t$ in a short interval $[t, t + \Delta t)$). We implement this Markov chain using the Gillespie algorithm [23] to obtain the simulated fixation probabilities in figures 5 and 6.</caption><thead><tr><th>event</th><th>jump<br>(ΔS<sub>u</sub><sup>n</sup>, ΔS<sub>v</sub><sup>n</sup>, ΔI<sub>uw</sub><sup>n</sup>, ΔI<sub>um</sub><sup>n</sup>, ΔI<sub>wv</sub><sup>n</sup>, ΔI<sub>vm</sub><sup>n</sup>, ΔR<sup>n</sup>)</th><th>rate</th></tr></thead><tbody><tr><td>birth</td><td>(1, 0, 0, 0, 0, 0)</td><td>nv</td></tr><tr><td>vaccination</td><td>(-1, 1, 0, 0, 0, 0)</td><td>vs<sub>u</sub><sup>n</sup></td></tr><tr><td>loss of immunity</td><td>(1, -1, 0, 0, 0, 0)</td><td>ω<sub>v</sub>S<sub>v</sub><sup>n</sup></td></tr><tr><td></td><td>(1, 0, 0, 0, 0, 0, -1)</td><td>ω<sub>v</sub>R<sup>n</sup></td></tr><tr><td rowspan="4">infection</td><td>(-1, 0, 1, 0, 0, 0, 0)</td><td>β<sub>w</sub> (I<sub>wv</sub><sup>n</sup> + I<sub>wv</sub><sup>n</sup>) S<sub>u</sub><sup>n</sup></td></tr><tr><td>(-1, 0, 0, 1, 0, 0, 0)</td><td>β<sub>m</sub> (I<sub>um</sub><sup>n</sup> + I<sub>vm</sub><sup>n</sup>) S<sub>u</sub><sup>n</sup></td></tr><tr><td>(0, -1, 0, 0, 1, 0, 0)</td><td>ε<sub>w</sub>β<sub>w</sub> (I<sub>wv</sub><sup>n</sup> + I<sub>wv</sub><sup>n</sup>) S<sub>v</sub><sup>n</sup></td></tr><tr><td>(0, -1, 0, 0, 0, 1, 0)</td><td>ε<sub>m</sub>β<sub>m</sub> (I<sub>um</sub><sup>n</sup> + I<sub>vm</sub><sup>n</sup>) S<sub>v</sub><sup>n</sup></td></tr><tr><td rowspan="4">recovery</td><td>(0, 0, -1, 0, 0, 0, 1)</td><td>γ<sub>w</sub>I<sub>wv</sub><sup>n</sup></td></tr><tr><td>(0, 0, 0, -1, 0, 0, 1)</td><td>γ<sub>m</sub>I<sub>um</sub><sup>n</sup></td></tr><tr><td>(0, 0, 0, 0, -1, 0, 1)</td><td>γ<sub>w</sub>I<sub>wv</sub><sup>n</sup></td></tr><tr><td>(0, 0, 0, 0, 0, -1, 1)</td><td>γ<sub>m</sub>I<sub>vm</sub><sup>n</sup></td></tr><tr><td rowspan="6">death</td><td>(-1, 0, 0, 0, 0, 0, 0)</td><td>δS<sub>u</sub><sup>n</sup></td></tr><tr><td>(0, -1, 0, 0, 0, 0, 0)</td><td>δS<sub>v</sub><sup>n</sup></td></tr><tr><td>(0, 0, -1, 0, 0, 0, 0)</td><td>(δ + α<sub>w</sub>)I<sub>wv</sub><sup>n</sup></td></tr><tr><td>(0, 0, 0, -1, 0, 0, 0)</td><td>(δ + α<sub>m</sub>)I<sub>um</sub><sup>n</sup></td></tr><tr><td>(0, 0, 0, 0, -1, 0, 0)</td><td>(δ + α<sub>w</sub>)I<sub>wv</sub><sup>n</sup></td></tr><tr><td>(0, 0, 0, 0, 0, -1, 0)</td><td>(δ + α<sub>m</sub>)I<sub>vm</sub><sup>n</sup></td></tr><tr><td></td><td>(0, 0, 0, 0, 0, 0, -1)</td><td>δR<sup>n</sup></td></tr></tbody></table>

viral adaptation are consistent with a broad range of acute infections of humans (e.g. SARS-CoV, influenza, measles; see table 2). In §4, we explore the robustness of our results after relaxing some of our simplifying assumptions.

## 2.1. Two approximations

Following [31], our analysis makes use of two approximations to our Markov process model. The first, deterministic approximation, uses ordinary differential equations (ODEs) and is appropriate when all types of host are abundant, but fails to correctly capture the dynamics when one or more types is rare (e.g. at the time of introduction of the mutant strain). The second uses a birth-and-death process (e.g. [32]) to approximate rare types and analyse stochastic phenomena, like extinction.

### 2.1.1. Deterministic approximation

For our first, deterministic approximation, we work with host densities defined by

$$X_i^n = S_i^n / n, \quad Y_{ij}^n = I_{ij}^n / n, \quad \text{and} \quad Z^n = R^n / n \qquad (2.2)$$

(i = u, v, j = w, m) and set

$$N^n = H^n / n = \sum_{i \in \{u, v\}} X_i^n + \sum_{\substack{i \in \{u, v\} \\ j \in \{w, m\}}} Y_{ij}^n + Z^n. \qquad (2.3)$$
For simplicity we assume that recovered hosts have lifelong immunity and thus reinfection is not possible. We assume frequency-dependent transmission where the number of contacts a host may have in the population is constant but a proportion of those contacts may be infectious. Note, however, that other forms of transmission (e.g. density-dependent transmission [33]) are expected to yield qualitatively similar results. We summarize the states of the process and the jump rates at which individuals transition between states in Table S1 and in Figure 2.

We use this model to examine the epidemiological and evolutionary dynamics following the start of the vaccination campaign. In the following, for the sake of simplicity, we focus on scenarios where we assume the pathogen population has reached an endemic equilibrium before the start of vaccination (but we consider alternative scenarios in the Discussion).

## 3 Results

### 3.1 Deterministic Approximation

As a first step in our analysis, we use a deterministic approximation for large values of $n$ [29] and we work with host densities defined as $S_i = S_i^n/n$, $I_{ij} = I_{ij}^n/n$, $N = N^n/n$. This corresponds to replacing discrete individuals by densities and interpreting the rates in Figure 2 as describing continuous flows rather than jumps. This yields a system of ordinary differential equations:

$$ 
\begin{aligned}
\dot{S}_U &= \lambda(1-p) + \omega S_V - \left( \beta_w \frac{I_{Uw} + I_{Vw}}{N} + \beta_m \frac{I_{Um} + I_{Vm}}{N} + \delta + \nu \right) S_U \\
\dot{S}_V &= \lambda p + \nu S_U - \left( e_w \beta_w \frac{I_{Uw} + I_{Vw}}{N} + e_m \beta_m \frac{I_{Um} + I_{Vm}}{N} + \delta + \omega \right) S_V \\
\dot{I}_{Uw} &= \beta_w (I_{Uw} + I_{Vw}) \frac{S_U}{N} - (\delta + \alpha_w + \gamma_w) I_{Uw} \\
\dot{I}_{Um} &= \beta_m (I_{Um} + I_{Vm}) \frac{S_U}{N} - (\delta + \alpha_m + \gamma_m) I_{Um} \\
\dot{I}_{Vw} &= e_w \beta_w (I_{Uw} + I_{Vw}) \frac{S_V}{N} - (\delta + \alpha_w + \gamma_w) I_{Vw} \\
\dot{I}_{Vm} &= e_m \beta_m (I_{Um} + I_{Vm}) \frac{S_V}{N} - (\delta + \alpha_m + \gamma_m) I_{Vm} \\
\dot{N} &= \lambda - \delta N - \alpha_w (I_{Uw} + I_{Vw}) - \alpha_m (I_{Um} + I_{Vm}),
\end{aligned}
\quad (1) $$
<table><caption>**Table 2.** Parameters and dynamical variables of the model.</caption><thead><tr><th>parameters<br>symbol</th><th>description</th><th>default value</th></tr></thead><tbody><tr><td>ν</td><td>influx rate of susceptible hosts</td><td>3 × 10<sup>-4</sup> week<sup>-1</sup></td></tr><tr><td>δ</td><td>natural death rate of hosts</td><td>3 × 10<sup>-4</sup> week<sup>-1</sup></td></tr><tr><td>ω<sub>r</sub></td><td>rate of waning immunity of naturally immune hosts</td><td>0.05 week<sup>-1</sup></td></tr><tr><td>ω<sub>v</sub></td><td>rate of waning immunity of vaccinated hosts</td><td>0.05 week<sup>-1</sup></td></tr><tr><td>υ</td><td>rate of vaccination</td><td>variable</td></tr><tr><td>Θ<sub>u</sub></td><td>rate of mutation (from w to m) in unvaccinated hosts</td><td>variable</td></tr><tr><td>Θ<sub>v</sub></td><td>rate of viral mutation (from w to m) in vaccinated hosts</td><td>variable</td></tr><tr><td>n</td><td>system size (scaling parameter allowing us to manipulate the pathogen population size)</td><td>variable</td></tr><tr><td>α<sub>w</sub></td><td>virulence (increased mortality rate) by strain w</td><td>0.02 week<sup>-1</sup></td></tr><tr><td>β<sub>w</sub></td><td>rate of transmission of strain w</td><td>10 week<sup>-1</sup></td></tr><tr><td>γ<sub>w</sub></td><td>recovery rate of the host infected by strain w</td><td>2 week<sup>-1</sup></td></tr><tr><td>ε<sub>w</sub></td><td>probability of infection of vaccinated hosts by strain w</td><td>0.05</td></tr><tr><td>R<sub>w</sub></td><td>reproduction number of strain w</td><td>4.95</td></tr><tr><td>α<sub>m</sub></td><td>virulence (increased mortality rate) by strain m</td><td>variable</td></tr><tr><td>β<sub>m</sub></td><td>rate of transmission of viral strain m</td><td>variable</td></tr><tr><td>γ<sub>m</sub></td><td>recovery rate of the host infected by strain m</td><td>2 week<sup>-1</sup></td></tr><tr><td>ε<sub>m</sub></td><td>probability of infection of vaccinated hosts by strain m</td><td>1</td></tr></tbody><thead><tr><th>variables<br>number</th><th>density</th><th></th></tr></thead><tbody><tr><td>S<sup>n</sup><sub>u</sub></td><td>unvaccinated susceptible hosts</td><td>X<sub>u</sub></td></tr><tr><td>S<sup>n</sup><sub>v</sub></td><td>vaccinated susceptible hosts</td><td>X<sub>v</sub></td></tr><tr><td>R<sup>n</sup></td><td>recovered hosts</td><td>Z</td></tr><tr><td>I<sup>n</sup><sub>w</sub></td><td>unvaccinated hosts infected by strain i ∈ {w, m}</td><td>γ<sub>w</sub></td></tr><tr><td>I<sup>n</sup><sub>v</sub></td><td>vaccinated hosts infected by strain i ∈ {w, m}</td><td>γ<sub>v</sub></td></tr></tbody></table>

As $n$ becomes large, the changes in the densities due to jumps in the Markov chain become smaller and smaller. As $n \to \infty$, the $X_{ij}^n$, $Y_{ij}^n$ and $Z^n$ approach limits $X_i$, $Y_{ij}$ and $Z$. These limits obey a system of ordinary differential equations:

$$
\begin{align}
\dot{X}_u &= \nu + \omega_v X_v + \omega_r Z - \left( \beta_w \frac{Y_{uw} + Y_{vw}}{N} + \beta_m \frac{Y_{um} + Y_{vm}}{N} + \delta + \nu \right) X_u, \nonumber \\
\dot{X}_v &= \nu X_u - \left( \epsilon_w \beta_w \frac{Y_{uw} + Y_{vw}}{N} + \epsilon_m \beta_m \frac{Y_{um} + Y_{vm}}{N} + \delta + \omega_v \right) X_v, \nonumber \\
\dot{Y}_{uw} &= \beta_w (Y_{uw} + Y_{vw}) \frac{X_u}{N} - (\delta + \alpha_w + \gamma_w) Y_{uw}, \nonumber \\
\dot{Y}_{um} &= \beta_m (Y_{um} + Y_{vm}) \frac{X_u}{N} - (\delta + \alpha_m + \gamma_m) Y_{um}, \tag{2.4} \\
\dot{Y}_{vw} &= \epsilon_w \beta_w (Y_{uw} + Y_{vw}) \frac{X_v}{N} - (\delta + \alpha_w + \gamma_w) Y_{vw}, \nonumber \\
\dot{Y}_{vm} &= \epsilon_m \beta_m (Y_{um} + Y_{vm}) \frac{X_v}{N} - (\delta + \alpha_m + \gamma_m) Y_{vm}, \nonumber \\
\dot{Z} &= (\gamma_w Y_{uw} + \gamma_w Y_{vw} + \gamma_m Y_{um} + \gamma_m Y_{vm}) - (\delta + \omega_r) Z. \nonumber
\end{align}
$$

This corresponds to replacing discrete individuals by continuous densities and interpreting the rates in figure 2 as describing continuous flows rather than jumps (see example B on p. 453 and theorem 11.2.1 on p. 456 in [33] for the details and proofs of this approximation; chapter 5 in [34] gives a readable summary with an epidemiological focus).
92

It is also convenient to track the dynamics of the total density of hosts infected with the same strain i, $I_i := I_{Ui} + I_{Vi}$, which yields:

94

$$ \dot{I}_i = \underbrace{\left( \left( \beta_i \frac{S_U}{N} + e_i \beta_i \frac{S_V}{N} \right) - (\delta + \alpha_i + \gamma_i) \right)}_{r_i = \text{growth rate of strain i}} I_i \quad (2) $$

The ability of the strains to grow is given by the sign of the growth rate $r_i$. Note that this growth rate depends on the four different traits of the pathogen: $\alpha_i, \beta_i, \gamma_i, e_i$. But this growth rate depends also on the densities $S_U$ and $S_V$, which vary with $t$, the time since the start of vaccination (i.e., vaccination starts at $t = 0$). The coefficient of selection $s_m$ on the mutant strain relative to the wild type is:

96

97

98

99

$$ s_m = r_m - r_w = (\beta_m - \beta_w) \frac{S_U}{N} + (e_m \beta_m - e_w \beta_w) \frac{S_V}{N} - (\alpha_m - \alpha_w + \gamma_m - \gamma_w) \quad (3) $$

In other words, both the genetics (the traits of strain i) and the environment (the epidemiological state of the host population) govern selection and strain dynamics.

101

## 3.2 Pathogen eradication and vaccination threshold

The ability of the wild-type pathogen population to grow can be measured by its effective per-generation reproduction ratio which is given by:

103

104

$$ R_w^e = R_w \left( \frac{S_U}{N} + e_i \frac{S_V}{N} \right) \quad (4) $$

105

106

107

108

109

110

where $R_i = \frac{\beta_i}{\delta+\alpha_i+\gamma_i}$. Hence, a reduction of the availability of susceptible hosts with vaccination may drive down the density of the wild-type pathogen when the production of new infected hosts (infection “birth”) does not compensate for the recovery and death of infected hosts (infection “death”), that is when $R_w^e < 1$. Ultimately, vaccination can even lead to the eradication of the wild-type pathogen (Figure 1a) either when the vaccine is sufficiently efficient ($e_w R_w > 1$) or when the vaccination coverage is sufficiently high [34, 16]. The deterministic model (1) can be used
It is also convenient to track the dynamics of the total density of hosts infected with the same strain i, $Y_i := Y_{ui} + Y_{vi}$, which yields

$$ \dot{Y}_i = \underbrace{\left( \left( \beta_i \frac{X_u}{N} + \epsilon_i \beta_i \frac{X_v}{N} \right) - (\delta + \alpha_i + \gamma_i) \right)}_{r_i = \text{growth rate of strain i}} Y_i. \quad (2.5) $$

The ability of the strain i to grow is given by the sign of the growth rate $r_i$. Note that this growth rate depends on the four different traits of the pathogen: $\alpha_i, \beta_i, \gamma_i, \epsilon_i$. The growth rate also depends on the densities $X_u(t)$ and $X_v(t)$, which vary with $t$, the time since the start of vaccination (i.e. vaccination starts at $t=0$). For simplicity, we assume that at time $t=0$, the wild-type is at its endemic equilibrium (see electronic supplementary material, SI.1, for details), and that there are no mutants (we relax this assumption in the electronic supplementary material, §5).

The coefficient of selection, $s_m(t)$, of the mutant strain relative to the wild-type is

$$ s_m(t) = r_m(t) - r_w(t) = (\beta_m - \beta_w) \frac{X_u(t)}{N(t)} + (\epsilon_m \beta_m - \epsilon_w \beta_w) \frac{X_v(t)}{N(t)} - (\alpha_m - \alpha_w + \gamma_m - \gamma_w). \quad (2.6) $$

In other words, both the genetics (the phenotypic traits of strain i) and the environment (the epidemiological state of the host population) govern selection and strain dynamics.

### 2.1.2. Birth-and-death process approximation

Suppose that a mutant strain appears at time $t_{int} \ge 0$ in a single infected host, $I_m^n(t_{int}) = 1$, that is, with density $Y_m^n(t_{int}) = \frac{1}{n}$. Taking $n \to \infty$, we get $Y_m(t_{int}) = 0$. Using this as an initial condition in (2.4), we find that $Y_m(t) \equiv 0$ for all $t \ge t_{int}$. This does not mean that the mutant is absent, but rather that it is not yet sufficiently abundant to be visible at the coarse resolution of the ODE approximation, (2.4). In particular, while rare, the mutant strain does not have a detectable effect on the density of susceptible hosts.

To account for the rare mutant, we use (2.4) to define a birth-and-death process, $\tilde{I}_m(t)$, that approximates the number of individuals infected with the mutant strain at times $t \ge t_{int}$, and allows us to estimate the probabilities of invasion (§3.2.2) and fixation (§5.1) of the mutant strain.

Each death in the birth-and-death process corresponds to the removal of an infected individual, which occurs by host death or recovery at combined rate

$$ d_m = \delta + \alpha_m + \gamma_m. \quad (2.7) $$

We approximate the rate of new infections,

$$ \frac{\beta_m(S_u^n(t) + \epsilon_m S_v^n(t))}{H^n(t)} = \frac{\beta_m(X_u^n(t) + \epsilon_m X_v^n(t))}{N^n(t)}, $$

by replacing the stochastic quantities $X_u^n(t)$, $X_v^n(t)$ and $N^n(t)$ by their deterministic approximations $X_u(t)$, $X_v(t)$ and $N(t)$, giving the time-dependent birth (i.e. 'birth' of an infection) rate:

$$ b_m(t) = \frac{\beta_m(X_u(t) + \epsilon_m X_v(t))}{N(t)}. \quad (2.8) $$

As we observed above, for the deterministic approximation, $Y_m(t_{int}) = 0$, and so we can compute $X_u(t)$, $X_v(t)$ and $N(t)$ using (2.4) without the mutant strain, using initial conditions (electronic supplementary material, SI.1). See [17, electronic supplementary material, §8.2] for a rigorous justification for this approximation.

The so-called 'merciless dichotomy' [35] tells us that, starting with one individual, the birth-and-death process either goes extinct, or grows indefinitely. Thus, either the mutant strain vanishes, or the number infected with the mutant strain will eventually grow to be of the order of $n$ individuals, after which we can use (2.4) to compute the densities of both wild-type and mutant strains. Arguing as above, we can compute the probability that either the wild-type or mutant strain goes extinct when susceptible hosts are depleted at the end of a major epidemic, using birth-and-death processes $\tilde{I}_w(t)$ and $\tilde{I}_m(t)$, provided we choose suitable initial conditions for (2.4) and the birth-and-death processes to reflect the conditions at the end of a major epidemic; see [31]. Indeed, provided we correctly choose initial conditions for [4] to account for successful invasion of the mutant—see electronic supplementary material, §4—we can use the birth-and-death processes throughout the epidemic phase: when infected hosts are abundant, the infection and recovery rates for both ODE and birth-and-death process approximations coincide, so the law of large numbers ensures that the densities (i.e. $\frac{\tilde{I}_w(t)}{n}$ and $\frac{\tilde{I}_m(t)}{n}$) predicted via the birth-and-death process approximation will agree with those from the ODE approximation when $n$ is large.

## 3. Results

### 3.1. Pathogen eradication and vaccination threshold

The ability of strain i to increase can be measured by its effective per-generation reproduction ratio:

$$ \mathcal{R}_i^e(t) = \mathcal{R}_i \left( \frac{X_u(t)}{N(t)} + \epsilon_i \frac{X_v(t)}{N(t)} \right), \quad (3.1) $$
to identify the threshold $\nu_c$ of the speed of vaccination rollout above which the wild-type pathogen can be driven to extinction (see Methods and SI Section 1):

$$ \nu_c = \frac{R_w(\delta(1 - (1 - e_w)p) + \omega) - (\delta + \omega)}{1 - R_w e_w} \quad (5) $$

As expected, better vaccines (i.e., lower values of $e_w$ and $\omega$) yield lower threshold values for the speed of vaccination. Imperfect vaccines (i.e., higher values of $e_w$ and $\omega$), in contrast, are unlikely to allow eradication. Note that, if we wait sufficiently long, the population of the wild-type pathogen will be driven to extinction in a stochastic way even when $\nu < \nu_c$. Indeed, in a finite host population, sooner or later, the pathogen population is doomed to become extinct because of demographic stochasticity, but the extinction time when $\nu < \nu_c$ will usually be very long, increasing exponentially with the system size $n$. From now on, we are going to neglect the possibility of extinction of the wild type due to vaccination when $\nu < \nu_c$ (which is a good approximation when $n$ is large) but will return to it in the discussion.

The spread of a new pathogen variant may erode the efficacy of vaccination and, consequently, could affect the ability to control and, ultimately, to eradicate the pathogen. But before the replacement of the wild type by a vaccine-escape variant the pathogen population may go through three steps that may result (or not) in pathogen adaptation: (1) introduction of the vaccine-escape variant by mutation, (2) extinction (Figure 1c) or invasion (Figure 1d-e) of the vaccine-escape variant, (3) fixation (Figure 1e) or not of the vaccine-adapted variant. Each of these steps is very sensitive to stochasticity because the number of vaccine-escape variants is very small in the early phase of its emergence.

### 3.3 Step 1: Introduction of the variant by mutation

The first step of adaptation is driven by the production of new variants by the wild-type pathogen through mutation. The level of adaptation to unvaccinated and vaccinated hosts may vary among those variants [10]. Vaccine-escape mutations that do not carry any fitness costs (or may even be adaptive) in unvaccinated hosts are expected to invade and fix relatively easily irrespective of
where $\mathcal{R}_i = \frac{\beta_i}{\delta + \alpha_i + \gamma_i}$, i = m, w. Hence, a reduction of the availability of susceptible hosts with vaccination may drive down the density of the wild-type pathogen when the production of new infected hosts (infection 'birth') does not compensate for the recovery or death of infected hosts (infection 'death'), that is, when $\mathcal{R}_w^e < 1$. Ultimately, vaccination can even lead to the eradication of the wild-type pathogen (figure 1a) either when the vaccine is sufficiently efficient ($\varepsilon_w \mathcal{R}_w < 1$) or when the vaccination coverage is sufficiently high [1,2].

Interestingly, if the aim is to eradicate an already established disease, bringing the reproduction number of the wild-type strain at the disease free equilibrium below one (i.e. $\mathcal{R}_w^\infty < 1$ (see electronic supplementary material, SI.3)) may not be sufficient to do so. Indeed, as pointed out by several earlier studies [36,37], imperfect vaccination may yield a backward bifurcation at the disease free equilibrium. In this case, the pathogen may persist even when vaccination brings $\mathcal{R}_w^\infty$ below one. Nonetheless, an analysis of our model indicates that the condition for the emergence of backward bifurcation are very limited (see electronic supplementary material, §1.3). In the following, we thus use the condition $\mathcal{R}_w^\infty < 1$ to identify the critical rate $v_c$ of the rate of vaccination above which the wild-type pathogen can be driven to extinction (see electronic supplementary material, §1.3):

$$v_c = \frac{(\mathcal{R}_w - 1)(\delta + \omega_v)}{1 - \mathcal{R}_w \varepsilon_w} \quad (3.2)$$

As expected, better vaccines (i.e. lower values of $\varepsilon_w$ and $\omega_v$) yield lower threshold values for the rate of vaccination. Imperfect vaccines (i.e. higher values of $\varepsilon_w$ and $\omega_v$), by contrast, are unlikely to allow eradication. Note that the population of the wild-type pathogen will be eventually be driven to extinction by stochastic fluctuations even when $v < v_c$ [38,39]. Indeed, in a finite host population, sooner or later, the pathogen population is doomed to go extinct because of demographic stochasticity, but the extinction time when $v < v_c$ will usually be very long, increasing exponentially with the system size $n$ [40–42]. From now on, we neglect the possibility of extinction of the wild-type due to vaccination when $v < v_c$ (which is a good approximation when $n$ is large).

The spread of a new pathogen variant may erode the efficacy of vaccination and, consequently, could affect the ability to control and, ultimately, to eradicate the pathogen. However, before the replacement of the wild-type by a vaccine-escape variant the pathogen population may go through three steps that may ultimately result (or not) in pathogen adaptation: (i) introduction of the vaccine-escape variant by mutation, (ii) extinction (figure 1c) or invasion (figure 1d–f) of the vaccine-escape variant introduced by mutation, (iii) fixation (figure 1f) or not (figure 1d,e) of the invading vaccine-escape variant. Each of these steps is very sensitive to demographic stochasticity as the number of vaccine-escape variants is very small in the early phase of its emergence. This motivates our birth-and-death process approximation (§2.1.2).

## 3.2. The steps of pathogen adaptation

Using the two approximations above (§2.1), we quantify the steps of pathogen evolution. First, we consider the appearance of a novel vaccine-resistant variant, which will either rapidly go extinct, or invade, causing an epidemic outbreak. Then, at the end of an epidemic, susceptible hosts are depleted, there are few individuals remaining infected with either wild-type or mutant strains, and both strains are at risk of extinction. If the variant outlives the wild-type, then the pathogen has adapted to the vaccine.

### 3.2.1. Step 1: introduction of the variant by mutation

The first step of adaptation is driven by the production of new variants of the wild-type pathogen through mutation. The degree of adaptation to unvaccinated and vaccinated hosts may vary among those variants [9]. For instance, some vaccine-escape mutations may carry no fitness costs (or may even be adaptive) in unvaccinated hosts. These variants would be expected to invade and fix because they are strongly favoured by natural selection when the proportion of vaccinated hosts builds up. They will have a high probability of avoiding the risk of early extinction irrespective of the vaccination strategy. We thus focus on variants that carry fitness costs in immunologically naive hosts (i.e. variants specialized on vaccinated hosts [9]). In principle, the introduction of the vaccine-escape mutation may occur before the rollout of vaccination. The distribution of these mutations is expected to follow a stationary distribution resulting from the action of recurrent mutations and negative selection (see electronic supplementary material, §5). If the fitness cost in naive hosts is high, and/or if the mutation rate is low, then these pre-existing mutants are expected to be rare. In the following, we neglect the presence of pre-existing mutants and we focus on a scenario where the first vaccine-escape mutant appears after the start of vaccination (but see electronic supplementary material, §5, where we discuss the effect of standing genetic variation).

At the onset of the vaccination campaign (i.e. $t = 0$) we assume that the system is at the endemic equilibrium (the equilibrium densities $X_u(0)$, $Y_{uw}(0)$ and $Y_{vw}(0)$ are given in electronic supplementary material, SI.1). We assume that an individual host infected with the wild-type produces vaccine-escape mutants at a small, constant rate $\mu_u$ if unvaccinated and $\mu_v$ if vaccinated. We make the standard population genetic assumption (e.g. [43, chapter 5]) that

$$0 < \theta_u = \lim_{n \to \infty} n \mu_u = \mathcal{O}(1) \quad \text{and} \quad 0 < \theta_v = \lim_{n \to \infty} n \mu_v = \mathcal{O}(1), \quad (3.3)$$

which ensures that the number of adaptive mutations remains non-zero and bounded in the large population limit that gives the deterministic approximation (2.4). In practice, we use the infinite-population limit to approximate finite populations, and we expect fewer adaptive mutations in a given period when modelling e.g. a single city versus an entire nation. We thus have to modify the values of $\theta_u$ and $\theta_v$ to reflect the size of the community. While the rate of mutation is assumed to be constant through time, whether or not a mutant will escape extinction within a host may depend on the type of host. Indeed, a vaccine-escape mutation
the vaccination strategy. We thus focus on variants that carry fitness costs in immunologically naïve hosts (i.e., variants specialized on vaccinated hosts [10]). In principle, the introduction of the vaccine-escape mutation may occur before the rollout of vaccination. The distribution of these mutations is expected to follow a stationary distribution resulting from the action of recurrent mutations and negative selection (see Methods). If these fitness costs are high and/or if the mutation rate is low these pre-existing mutants are expected to be rare. In the following, we neglect the presence or pre-existing mutants and we focus on a scenario where the first vaccine-escape mutant appears after the start of vaccination (but see Methods, section 5.5).

At the onset of the vaccination campaign (i.e., $t = 0$) we assume that the system is at the endemic equilibrium (i.e., the equilibrium densities $S_U^0$, $I_{Uw}^0$ and $I_{Vw}^0$ are given in (11) in the Methods). We assume that an individual host infected with the wild type produces vaccine-escape mutants at a small, constant rate $\theta_U/n$ if unvaccinated and $\theta_V/n$ if vaccinated. In addition, we assume that $\theta_U$ and $\theta_V$ are small enough that within-host clonal interference is unlikely, and that $\theta_V \ge \theta_U$ to account for the fact that within-host selection may favor mutants in vaccinated hosts. The total production of mutants is thus equal to $\frac{\theta_U}{n}I_{Uw}^n(t) + \frac{\theta_V}{n}I_{Vw}^n(t) \approx \theta_U I_{Uw}(t) + \theta_V I_{Vw}(t)$, so that the probability density $f_m(t)$ of the arrival time $t = t_m$ of the vaccine-escape mutant is approximated by:

$$f_m(t) = (\theta_U I_{Uw}(t) + \theta_V I_{Vw}(t)) e^{-\int_0^t (\theta_U I_{Uw}(s) + \theta_V I_{Vw}(s)) ds}. \quad (6)$$

In other words, the time $t_m$ at which the vaccine-escape variant is introduced by mutation depends on the dynamics of the incidence of the infections by the wild type. Plots of $f_m$ for different values of rollout speed $\nu$ in **Figure 3** show that a faster rollout of vaccination delays the introduction of the vaccine-escape mutant. Indeed, a faster rollout is known to result in a drop of the incidence (the honey-moon period) [34, 13] which is expected to decrease the influx of new variants during this period. **Figure 3** also shows how higher values of $\theta_V$ can increase the influx of vaccine-escape variants. As discussed in the following section, the subsequent fate of vaccine-escape mutants depends strongly on the timing of their arrival.
may have a higher probability to escape within-host extinction in vaccinated hosts. We account for this effect by making a distinction between $\mu_u$ and $\mu_v$. If vaccine-escape mutations are more likely to escape extinction in vaccinated hosts we expect $\mu_v > \mu_u$. In other words, $\mu_v/\mu_u - 1$ is a measure of the within-host fitness advantage of the vaccine-escape mutant in vaccinated hosts (they are assumed to have the same within-host fitness in naïve hosts). We assume that $\mu_u$ and $\mu_v$ are small enough that within-host clonal interference among vaccinated-adapted variants is negligible. The total rate of production of mutants is thus equal to

$$ \mu_u I_{uw}^n(t) + \mu_v I_{vw}^n(t) \approx \theta_u Y_{uw}(t) + \theta_v Y_{vw}(t). \quad (3.4) $$

The arrival times of novel mutants are thus well approximated by a non-homogeneous Poisson process [44, p. 4] with rate

$$ \lambda_{\text{int}}(t) = \theta_u Y_{uw}(t) + \theta_v Y_{vw}(t). \quad (3.5) $$

The cumulative density function of the arrival time $T_{\text{int}}$ of the first vaccine-escape mutant is thus approximated by

$$ F_{\text{int}}(t) = \mathbb{P}\{T_{\text{int}} \le t\} = 1 - e^{-\int_{0}^{t} \lambda_{\text{int}}(s) ds}. \quad (3.6) $$

In other words, the time $T_{\text{int}}$ at which the vaccine-escape variant is first introduced by mutation depends on the dynamics of the incidence of the infections by the wild-type. Plots of $F_{\text{int}}(t)$ for different values of rollout speed $v$ in figure 3 show that a faster rollout of vaccination delays the introduction of the vaccine-escape mutant. This effect is particularly marked when $\omega_r = 0$ because life-long immunity is known to result in a massive transient drop of the incidence (the honeymoon period) [2,45] which is expected to decrease the influx of new variants during this period (electronic supplementary material, figure SI.1). Figure 3 also shows how higher values of $\omega_v$ can increase the influx of vaccine-escape variants. As discussed in the following section, the subsequent fate of vaccine-escape mutants depends strongly on the timing of their arrival.

### 3.2.2. Step 2: variant invasion

Immediately after its introduction, the dynamics of the vaccine-escape mutant may be approximated by a time-inhomogeneous birth-and-death process where the rate of birth (i.e. rate of new infections by the mutant) varies with the availability of susceptible hosts (see §2.1.2). The probability a mutant introduced at $T_{\text{int}} = t_{\text{int}}^1$ successfully invades (see [46] and electronic supplementary material, §2) is

$$ P_{\text{inv}}(t_{\text{int}}) = \frac{1}{1 + \int_{t_{\text{int}}}^{\infty} d_m e^{-\int_{t_{\text{int}}}^{s} b_m(u) - d_m du} ds}, \quad (3.7) $$

with $b_m(t)$ and $d_m$ as defined above, (2.7), (2.8). In general, the integrals in (3.7) are impossible to compute exactly. In §5.2.2, we describe a fast numerical method that allows us to compute $P_{\text{inv}}(t_{\text{int}})$.

Plotting the probability of invasion against the time of introduction, $t$, in figure 4 shows that the time at which the vaccine-escape mutant is introduced has a dramatic impact on the probability of escaping early extinction. If the mutant is introduced early, the density of susceptible vaccinated hosts remains very low, and the selection for the vaccine-escape mutant is too small to prevent stochastic extinctions. The probability of invasion increases with selection, and thus with the density of vaccinated hosts, which tends to increase with time (see equation (2.6)).

Taking $t \to \infty$ allows us to consider the situation when the vaccine-escape mutant appears at the post-vaccination endemic equilibrium, i.e. when the densities of unvaccinated and vaccinated susceptible hosts are $X_u^*$ and $X_v^*$, respectively (see electronic supplementary material, §1.3). At that point in time the effective per-generation reproduction ratio of genotype $i$ (i.e. the expected number of secondary infections produced by pathogen genotype $i$) is (cf. (3.1))

$$ \mathcal{R}_i^* = \lim_{t \to \infty} \mathcal{R}_i^e(t) = \mathcal{R}_i \left( \frac{X_u^*}{N^*} + \epsilon_i \frac{X_v^*}{N^*} \right). \quad (3.8) $$

By definition, at the endemic equilibrium set by the wild-type pathogen we have $\mathcal{R}_w^* = 1$. Hence, a necessary condition for the mutant to invade this equilibrium is $\mathcal{R}_m^* > 1$, i.e. the effective reproduction number of the mutant must be higher than that of the wild-type (see electronic supplementary material, §1.3). However, this is not a sufficient condition: many mutants that satisfy this condition will nonetheless rapidly go extinct due to demographic stochasticity. But by contrast to an early introduction of the mutant discussed above, the stochastic dynamics of the mutant is approximately a time-homogeneous birth-and-death process because the birth/infection rate of the mutant approaches $b_m^* = \beta_m \left( \frac{X_u^*}{N^*} + \epsilon_m \frac{X_v^*}{N^*} \right)$. This rate is constant because the density of susceptible hosts remains constant at the endemic equilibrium. The probability of mutant invasion after introducing a single host infected by the mutant is thus (see electronic supplementary material §3; figure 4)

$$ P_{\text{inv}}^* = \lim_{t \to \infty} P_{\text{inv}}(t) = 1 - \frac{\mathcal{R}_w^*}{\mathcal{R}_m^*} = 1 - \frac{1}{\mathcal{R}_m^*} \quad (3.9) $$

(note that we recover the strong-selection result of [17]). This expression shows that *at this endemic equilibrium* the fate of the mutant is fully governed by the per-generation reproduction ratio of the two strains, but does not depend on the specific values of the life-history traits of the mutant (provided the different vaccine-escape variants have the same value of $\mathcal{R}_m^*$).

<footnote n="1">¹To clarify, $T_{\text{int}}$ is the random time at which a mutation arises; when we specify $T_{\text{int}} = t_{\text{int}}$, we are conditioning on the event in which the random quantity $T_{\text{int}}$ takes the fixed value $t_{\text{int}}$.</footnote>
## 3.4 Step 2: Variant invasion

Immediately after its introduction, the dynamics of the vaccine-escape mutant may be approximated by a time-inhomogeneous birth-death process where the rate of birth (i.e., rate of new infections by the mutant) varies with the availability of susceptible hosts (see Methods, section 5.3). The probability that a single mutant with time-varying birth rate $b_m(t) = \beta_m \left( \frac{S_U(t)}{N(t)} + \frac{e_m S_V(t)}{N(t)} \right)$ and constant death rate $d_m = \delta + \alpha_m + \gamma_m$, introduced at time $t_m$, successfully invades (see [24] and SI, Section 2) is:

$$P_m^{t_m} = \frac{1}{1 + \int_{t_m}^{\infty} d_m e^{-\int_{t_m}^{t} b_m(s) - d_m ds}} dt \quad (7)$$

Plotting the probability of invasion vs. $t_m$ in Figure 4 shows that the time at which the vaccine-escape mutant is introduced has a dramatic impact on the probability of escaping early extinction. If the mutant is introduced early, the density $S_V(t)$ of susceptible vaccinated hosts remains very low and the selection for the vaccine-escape mutant is too small to prevent stochastic extinctions. The probability of invasion increases with selection, and thus with the density of vaccinated hosts, which tends to increase with time (see equation (3)).

Taking $t_m \to \infty$ allows us to tackle the situation when the vaccine-escape mutant appears at the endemic equilibrium, i.e., when the densities of unvaccinated and vaccinated susceptible hosts are $S_U^*$ and $S_V^*$, respectively. At that point in time the effective per-generation reproduction ratio of genotype $i$ (i.e. the expected number of secondary infections produced by pathogen genotype $i$) is:

$$R_i^* = R_i \left( \frac{S_U^*}{N^*} + e_i \frac{S_V^*}{N^*} \right) \quad (8)$$

By definition, at the endemic equilibrium set by the wild-type pathogen we have $R_w^* = R_w^e = 1$. Hence, a necessary condition for the mutant to invade this equilibrium is $R_m^* > 1$, i.e., the effective reproduction number of the mutant has to be higher than that of the wild type (see SI, Section 1). However, this is not a sufficient condition: many mutants that satisfy this condition will rapidly go extinct due to demographic stochasticity. But in contrast to an early introduction of the mutant discussed above, the stochastic dynamics of the mutant is approximately a time-homogeneous branching process because the birth rate of the mutant approaches $b_m^* = \beta_i \left( \frac{S_U^*}{N^*} + e_i \frac{S_V^*}{N^*} \right)$. This
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Faster vaccine rollout delays the arrival time of the first escape mutant. We plot the probability, $F_{\text{int}}(t)$, that the first escape mutant arrives prior to time $t$ for different speeds of vaccination rollout: $v = 0.05$ (a), 0.15 (b) and 0.24 (c). We contrast a scenario where $\theta_v = \theta_u$ (dashed line), and $\theta_v = 10 \times \theta_u$ (full line). Other parameter values: $\theta_u = 1$, $v = \delta = 3 \times 10^{-4}$, $\omega_v = \omega_r = 0.05$, $\alpha_w = 0.02$, $\beta_w = 10$, $\gamma_w = 2$, $\epsilon_w = 0.05$, $\mathcal{R}_w = 4.95$. For these parameter values, the critical rate of vaccination $v_c$ above which the wild-type pathogen is driven to extinction is $v_c \approx 0.264$ (see equation (3.2)).</figcaption></figure>

Interestingly, unlike $P_{\text{inv}}^{*}$, the probability, $P_{\text{inv}}(t_{\text{int}})$, that a mutant introduced at time $T_{\text{int}} = t_{\text{int}}$ successfully invades (3.7) is not governed solely by $\mathcal{R}_{i}$, but rather depends on the life-history traits of the mutants. For instance, assume that two vaccine-escape mutants have the same values of $\mathcal{R}_{m}$ and $\epsilon_{m}$ but they have very different life-history strategies. The ‘slow’ strain has low rates of transmission and virulence (in green in figure 4) while the ‘fast’ strain has high rates of transmission and virulence (in red in figure 4). Figure 4 shows that the high mortality rate of hosts infected by the fast strain increases the risk of early extinction and lowers the probability of invasion relative to the slow strain. Hence, in the early stage of adaptation, pathogen life-history matters and favours slow strains with lower rates of transmission and virulence.
birth rate is constant because the density of susceptible hosts remains constant at the endemic equilibrium. The probability of mutant invasion after introducing a single host infected by the mutant is thus (see SI Section 1; Figure 4):

$$P_m^* = \lim_{t_m \to \infty} P_m^{t_m} = 1 - \frac{R_w^*}{R_m^*} = 1 - \frac{1}{R_m^*} \quad (9)$$

Note that we recover the strong-selection result of [36]. This expression shows that *at this endemic equilibrium* the fate of the mutant is fully governed by the per-generation reproduction ratio of the two strains, but does not depend on the specific values of the life-history traits of the mutant (provided the different vaccine-escape variants have the same value of $R_m^*$).

Interestingly, unlike $P_m^*$, the probability $P_m^{t_m}$ of mutant invasion at time $t_m$ given in (7) is not governed solely by $R_i$, but rather depends on the life-history traits of the mutants. For instance, assume that two vaccine-escape mutants have the same values of $R_m$ and $e_m$ but they have very different life-history strategies. The “slow” strain has low rates of transmission and virulence (in green in Figure 4) while the “fast” strain has high rates of transmission and virulence (in red in Figure 4). Figure 4 shows that the high mortality rate of hosts infected by the fast strain increases the risk of early extinction and lowers the probability of invasion relative to the slow strain. Hence, in the early stage of adaptation, pathogen life-history matters and favours slow strains with lower rates of transmission and virulence.

## 3.5 Step 3: After variant invasion

Successful invasion of the vaccine-escape mutant means that it escaped the “danger zone” when its density is so low that it is very likely to go extinct (Figure 1d-f). After this invasion we can describe the dynamics of the polymorphic pathogen population using the deterministic approximation (1).

Because the invasion of the mutant at the endemic equilibrium set by the wild type requires that $R_m^* > R_w^*$, we may expect from the analysis of the deterministic model that the mutant would always replace the wild-type pathogen. That is, the wild-type pathogen would go extinct before the mutant (Figure 1f). This is indeed the case when the phenotypes of the mutant and the wild type are not very different because of the “invasion implies fixation” principle [17, 5, 38]. Yet, this
Probability of invasion of vaccine-escape mutant

<figure><img src="image_9.png" /><figcaption>Arrival time $t_{int}$ of the vaccine-escape mutant</figcaption></figure>

**Figure 4.** Probability of invasion of the vaccine-escape mutant increases with $T_{int}$. We plot the probability invasion $P_{inv}(t_{int})$ of a *slow* (blue) and a *fast* (red) vaccine-escape mutant for different speeds of vaccination rollout: $v = 0.05$ (a), 0.15 (b) and 0.24 (c). The *slow* mutant: $\alpha_m = 0.02$, $\beta_m = 7$, $\gamma_m = 2$, $\varepsilon_m = 1$, $\mathcal{R}_m = 3.46$. The *fast* mutant: $\alpha_m = 4.0606$, $\beta_m = 21$, $\gamma_m = 2$, $\varepsilon_m = 1$, $\mathcal{R}_m = 3.46$. The probability of invasion $P_{inv}^*$ in the limit $t_{int} \to \infty$ (see equation (3.9)) is indicated with the dashed black line. Other parameter values as in figure 3: $v = \delta = 3 \times 10^{-4}$, $\omega_v = \omega_r = 0.05$, $\alpha_w = 0.02$, $\beta_w = 10$, $\gamma_w = 2$, $\varepsilon_w = 0.05$, $\mathcal{R}_w = 4.95$.

### 3.2.3. Step 3: after variant invasion

Successful invasion of the vaccine-escape mutant means that it escaped the 'danger zone' when its density is so low that it is very likely to go extinct (figure 1d-f). After this invasion, we can describe the dynamics of the polymorphic pathogen population using the deterministic approximation (2.4).

Because invasion of the mutant at the endemic equilibrium set by the wild-type requires that $\mathcal{R}_m^* > \mathcal{R}_w^*$, we might expect from the analysis of the deterministic model that the mutant would always replace the wild-type pathogen. That is, the wild-type pathogen will always go extinct before the mutant (figure 1f). This is indeed the case when the phenotypes of the mutant and the wild-type are not very different because of the 'invasion implies fixation' principle [47-49]. This principle may be violated, however, if the phenotype of the vaccine-escape mutant is very different from the phenotype of the wild-type.
211
212

principle may be violated if the phenotype of the vaccine-escape mutant is very different than the phenotype of the wild type.

213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238

First, if the two genotypes (*w* and *m*) are sufficiently specialized on the two hosts, the long-term coexistence of the two genotypes is possible (**Figure 1e**). Second, the vaccine-escape mutant may be driven to extinction before the wild type if its life-history traits induce massive epidemiological perturbations after its successful invasion (**Figure 1d**). As pointed out by previous studies, more transmissible and aggressive pathogen strategies may yield larger epidemics because the speed of the epidemic is governed by the per-capita growth rate *r*ᵢ, not the per-generation reproduction ratio *R*ᵢ [13]. This explosive dynamics is followed by a decline which results in a very low incidence of the vaccine-escape mutant. In a finite host population, this may result in the extinction of the vaccine-escape mutant before the wild type [42]. We capture this outcome with a hybrid analytical-numerical approach that computes the probability *P*<sub>fix</sub><sup>*t*<sub>m</sub></sup> that the wild type will go extinct before the mutant (see Methods, section 5.4). **Figure 5** shows that two vaccine-escape mutants may have very different probabilities of fixation, even if they have the same per-generation reproduction ratio. The numerical computation of the probability of fixation agrees very well with individual-based stochastic simulations (**Figure S1**). The faster strain is unlikely to go to fixation because invasion is followed by a period where the birth rate drops to very low levels (far below the mortality rates, **Figure S2**). In other words, a more aggressive strategy will more rapidly degrade its environment, by depleting susceptible hosts, which is known to increase the probability of extinction [6]. Interestingly, this effect is only apparent when the time of introduction *t*<sub>m</sub> is large. Indeed, when the mutant is introduced soon after the start of vaccination, its probability of invasion is already very low because its initial growth rate is negative (**Figure S2a, b, c**). When the mutant is introduced at intermediate times, the initial growth rate of the mutant is positive because some hosts are vaccinated (**Figure S2d, e, f**). If the vaccine-escape mutant is introduced later, the growth rate of the mutant is initially very high as many hosts are vaccinated (and thus susceptible to the vaccine-escape mutant) but this is rapidly followed by a drop in host density (especially pronounced with the faster strain) which prevents the long-term establishment of the faster strain (see **Figure S2g, h, i**).
First, the long-term coexistence of the two genotypes is possible (figure 1e). The coexistence requires that each genotype is specialized on distinct types of host. The wild-type is specialized on unvaccinate hosts (i.e. $R_w > R_m$) and the mutant is specialized on the vaccinated hosts (i.e. $\epsilon_m > \epsilon_w$). Intermediate rates of vaccination maintain a mix of vaccinated and unvaccinated host which promotes coexistence between the two genotypes (electronic supplementary material, figure SI.2). Second, the vaccine-escape mutant may be driven to extinction before the wild-type if its life-history traits induce large epidemiological perturbations after its successful invasion (figure 1d). As pointed out by previous studies, more transmissible and aggressive pathogen strategies may yield larger epidemics because the speed of the epidemic is governed by the per capita growth rate $r_i$, not by the per-generation reproduction ratio $R_i$ [45]. These explosive dynamics are driven by an over-exploitation of the host population, and are immediately followed by a large decline in the incidence of the vaccine-escape mutant. In a finite host population, this may result in the extinction of the vaccine-escape mutant before the wild-type [6]. We capture this outcome with a hybrid analytical-numerical approach that computes the probability, $P_{fix}(t_{int})$, that the wild-type will go extinct before a mutant introduced at time $T_{int} = t_{int}$ (see §5.1). The numerical computation of the probability of fixation agrees very well with individual-based stochastic simulations. figure 5 shows that two vaccine-escape mutants may have very different probabilities of fixation, even if they have the same per-generation reproduction ratio: a 'faster' strain, with high rate of transmission and elevated virulence, and a 'slower' strain, with reduced transmission rate, and lower virulence. The faster strain is unlikely to go to fixation because invasion is followed by a period where the birth rate drops to very low levels (far below the mortality rates; electronic supplementary material, figure SI.3). In other words, a more aggressive strategy will more rapidly degrade its environment, by depleting susceptible hosts, which is known to increase the probability of extinction [50]. Interestingly, this effect is only apparent when the time of introduction, $T_{int}$, is large. Indeed, when the mutant is introduced soon after the start of vaccination, its probability of invasion is already very low because its initial growth rate is negative (electronic supplementary material, figure SI.3a-c). When the mutant is introduced at intermediate times, the initial growth rate of the mutant is positive because some hosts are vaccinated (electronic supplementary material, figure SI.3d-f). If the vaccine-escape mutant is introduced later, the growth rate of the mutant is initially very high as many hosts are vaccinated (and thus susceptible to the vaccine-escape mutant) but this is rapidly followed by a drop in host density (especially pronounced with the faster strain) which prevents the long-term establishment of the faster strain (see electronic supplementary material, figure SI.3g-i).

### 3.2.4. The overall risk of pathogen adaptation

The overall probability that the pathogen will adapt to vaccination (i.e. that a vaccine-escape variant invades and eventually replaces or coexists with the wild-type) depends upon the probability that the mutation will arise (step 1) and the probability that this mutation will escape early extinction (step 2) and eventually go to fixation (step 3). It is particularly relevant to explore the effect of the rate of vaccination rollout on the overall probability that some vaccine-escape variant successfully invades at some time $T_{inv} \le t$ after the start of the vaccination campaign (steps 1 and 2, figure 6). Note that several variants can arise and fail to invade before finally some variant manages to invade. We can use the probability of invasion, $P_{inv}(t)$, of a variant introduced at time $t$ to characterize the distribution, $F_{inv}(t)$, of the first time, $T_{inv}$, at which a mutant is introduced that successfully invades. Using (3.5) and (3.9), this is

$$F_{inv}(t) = \mathbb{P}\{T_{inv} \le t\} = 1 - e^{-\int_{0}^{t} \lambda_{int}(s)P_{inv}(s) ds}. \quad (3.10)$$

Note that the probability that no vaccine-escape mutant will ever arise is

$$\mathbb{P}\{T_{int} = \infty\} = e^{-\int_{0}^{\infty} \lambda_{int}(s) ds}.$$

By contrast, the probability that no vaccine-escape mutant will ever invade is the larger probability

$$\mathbb{P}\{T_{inv} = \infty\} = e^{-\int_{0}^{\infty} \lambda_{int}(s)P_{inv}(s) ds}.$$

Moreover, $P_{inv}(t)$ converges as $t \to \infty$ to $P_{inv}^* = 1 - 1/R_m^*$, which is nonzero, so $\mathbb{P}\{T_{int} = \infty\} = 0$ if and only if

$$\int_{0}^{\infty} \lambda_{int}(s) ds = \infty.$$

This in turn is true if and only if $\mathbb{P}\{T_{inv} = \infty\} = 0$. That is, the probability of adaptation is 1 if and only if $\lambda_{int}(t)$ is not integrable. In other words, the probability of adaptation is 1 in the limit $t \to \infty$ when the wild-type is not driven to extinction by vaccination (i.e. $v < v_c$) which implies that there is an uninterrupted flux of mutation producing vaccine-escape variants. One of these mutants must eventually escape extinction and invade. Yet, the time needed for a successful variant to appear may be very long (equation (3.10) and figure 6).

When $v > v_c$, vaccination is expected to eradicate the disease rapidly in our model (but see electronic supplementary material, §1.3). However, an escape mutation may appear by mutation before eradication and rescue the pathogen population. This scenario fits squarely within the framework of classical 'evolutionary rescue' modelling [51-53]. Yet, vaccination rollout is unlikely to be fast enough to eradicate the wild-type pathogen and, in this case, the probability of adaptation goes to 1 when $t \to \infty$. Indeed, when $v < v_c$, a vaccine-escape variant will eventually appear by mutation and invade. But what is less clear is how fast this adaptation will take place. We can use equation (3.10) to explore the effect of the speed of adaptation on the probability of pathogen adaptation at time $t$ after the start of vaccination (i.e. the speed of adaptation). Crucially, the speed of pathogen adaptation is maximized for intermediate values of the speed of vaccination rollout. This is due to the antagonistic consequences the speed of the rollout
## 3.6 The overall risk of pathogen adaptation

The overall probability that the pathogen will adapt to vaccination (i.e. that a vaccine-escape variant will replace the wild-type) depends upon the probability that the mutation will arise (step 1) and the probability that this mutation will escape early extinction (step 2) and eventually go to fixation (step 3). It is particularly relevant to explore the effect of the speed of vaccination rollout on the overall probability that some vaccine-adapted variant invades before a time t after the start of the vaccination campaign (steps 1 and 2, Figure 6):

$$A_m^t = 1 - e^{-\int_0^t (\theta_U I_{Uw}(s) + \theta_V I_{Vw}(s)) P_m^s ds} \quad (10)$$

As expected, when $\nu < \nu_c$ this probability of adaptation goes to 1 when $t \to \infty$. Indeed, when vaccination cannot eradicate the wild type, a vaccine-adapted variant will eventually appear by mutation and invade. When $\nu > \nu_c$, we recover a classical evolutionary rescue scenario where the arrival and the spread of a vaccine-adapted variant may stop the eradication [32, 2, 4]. Since vaccination is unlikely to lead to eradication we focus here on a scenario where $\nu < \nu_c$ and we use equation (10) to analyse the effect of the speed of adaptation on the probability of pathogen adaptation at a time t after the start of vaccination. Crucially, the risk of pathogen adaptation is maximized for intermediate values of the speed of vaccination rollout. This is due to the antagonistic consequences the speed of the rollout has upon these two steps of adaptation (compare Figures 3 and 4). Faster rollout reduces the influx of new mutations, but increases selection for vaccine-escape mutations. Notice how the risk of adaptation drops with the speed of the rollout of vaccination before the critical speed, $\nu < \nu_c$, that may ultimately lead to eradication. When the speed of vaccination is just below $\nu_c$, vaccination coverage is not high enough to allow eradication, but it is high enough to reduce massively the probability of adaptation through the reduction of the influx of new mutations. As expected, variants with different life-history may have different probabilities of adaptation. Indeed, as pointed out in the previous section, faster variants have a lower probability of invasion when the variant is introduced at the beginning of the vaccination campaign (Figure 4). This effect is relatively small but it is expected to be magnified by the subsequent risk of
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Probability of fixation of the vaccine-escape mutant may be low when $T_{int}$ is large. We plot the probability of fixation of (a) a *slow* (blue) and (b) a *fast* (red) vaccine-escape mutant for an intermediate speed of vaccination rollout: $v = 0.15$. The *slow* mutant: $\alpha_m = 0.02, \beta_m = 7, \gamma_m = 2, \epsilon_m = 1, \mathcal{R}_m = 3.46$. The *fast* mutant: $\alpha_m = 4.0606, \beta_m = 21, \gamma_m = 2, \epsilon_m = 1, \mathcal{R}_m = 3.46$. The full coloured lines give the probability of fixation $P_{fx}(t_{inv})$ computed numerically (see §5.2.3) and the dots give the results of individual-based simulations (see §5.3) for different values of $n$ which affect the pathogen population size and the intensity of demographic stochasticity. We plot the probability of invasion $P_{inv}(t)$ (see figure 4) with dashed coloured line and its asymptotic value $P_{inv}^*$ with a dotted black line. Other parameter values as in figure 3: $v = \delta = 3 \times 10^{-4}, \omega_v = \omega_t = 0.05, p = 0, \alpha_w = 0.02, \beta_w = 10, \gamma_w = 2, \epsilon_w = 0.05, \mathcal{R}_w = 4.95$.</figcaption></figure>

has upon these two steps of adaptation (compare figures 5 and 6). Faster rollout reduces $\lambda_{int}$, the influx of new mutations, but increases $P_{inv}$ because higher vaccination coverage yield stronger selection for vaccine-escape mutations. Figure 6 illustrates how the speed of adaptation given in [18] between the time-varying probability, $F_{int}(t)$, that a variant is introduced by mutation before time $t$ and the probability, $P_{inv}(t)$, that this variant successfully invades. After an initial transient, $P_{inv}^*$ is a good approximation of this probability of invasion (see (3.9) and figures 6 and 2). We illustrate how $P_{inv}^*$ (and thus the risk of pathogen adaptation) varies with the rate of vaccination in (2.6).

## 4. Discussion

Vaccination is a powerful tool to control the spread of infectious diseases, but some pathogens evolve to escape the immunity triggered by vaccines (e.g. influenza, SARS-CoV-2). Will pathogens continue to adapt to the different vaccines that are being used to halt their spread? Does the likelihood of this adaptation depend on the speed of the vaccination rollout? To answer these questions we must first understand the different steps that may eventually lead to adaptation to vaccination.

Mutation is the fuel of evolution, and the first step of adaptation to vaccination is the mutational process that produces vaccine-escape variants. For instance, even if initial estimates of SARS-CoV-2 mutation rates were reassuringly low [54], the virus has managed to evolve higher rates of transmission [55,56] and these adaptations are challenging control measures currently being used to slow down the ongoing pandemic. The ability of the new variants of SARS-CoV-2 to escape immunity is also worrying and indicates that viral adaption can weaken vaccine efficacy [57,58]. The rate at which these potential vaccine-escape mutations are introduced depends on the density of hosts infected by the wild-type virus. In this respect, a faster rollout of vaccination is expected to delay the arrival of these mutations (figure 3). Some authors, however, have argued that vaccine-escape mutations may
extinction following invasion (Figure 5).

4 Discussion

Vaccination is a powerful tool to control the spread of infectious diseases, but some pathogens evolve to escape the immunity triggered by vaccines. Will SARS-CoV-2 continue to adapt to the different vaccines that are currently being used to halt the ongoing pandemic? Does the likelihood of this adaptation depend on the speed of the vaccination rollout? To answer these questions we must first understand the different steps that may eventually lead to adaptation to vaccination.

Mutation is the fuel of evolution, and the first step of adaptation to vaccination is the mutational process that produces vaccine-escape variants. Even if initial estimates of SARS-CoV-2 mutation rates were reassuringly low [39], the virus has managed to evolve higher rates of transmission [9, 46] and these adaptations are challenging current control measures used to slow down the ongoing pandemic. The ability of the new variants of SARS-CoV-2 to escape immunity is also worrying and indicates that viral adaption can weaken vaccine efficacy [47, 37]. The rate at which these potential vaccine-escape mutations are introduced depends on the density of hosts infected by the wild-type virus. In this respect, a faster rollout of vaccination is expected to delay the arrival of these mutations (Figure 3).

Some authors have argued that the emergence of vaccine-escape mutations may be more likely in infected hosts which are partially immunized [12, 8, 10]. Our model can be used to explore the consequences of this within-host evolution in vaccinated hosts (e.g., taking $θ_V > θ_U$). A larger value of $θ_V$ increases the overall rate of mutation (Figure 3) but this effect is modulated by the fraction of the host population that is vaccinated. Consequently, when $θ_V > θ_U$, the speed of vaccination rollout can have a non-monotonic effect on the probability that a vaccine-escape mutation is introduced (see Figure S3). Indeed, when the rate of vaccination remains low, the enhancing effect of vaccination on the rate of mutation can counteract the delaying effect of faster vaccination rollout discussed above. But the probability that a vaccine-escape mutation is introduced drops to very low levels when the rate of vaccination becomes overwhelmingly high.

The second step of adaptation starts as soon as the vaccine-escape mutant has been introduced
<figure><img src="image_11.png" /><figcaption>**Figure 6.** The probability of adaptation is maximized for intermediate speed of vaccination rollout. In (a), we plot the probability of adaptation $F_{\text{inv}}(t)$ (black lines) against the speed of vaccination rollout at different points in time. In (b), we plot the probability $F_{\text{int}}(t)$ of the introduction of at least one mutant before different points in time $t$ (orange lines) and the probability $P_{\text{inv}}^*(\text{purple line})$ which gives a good approximation of the probability of successful invasion of an escape-mutant. The dashed purple line gives the probability of invasion of the escape-mutant in the absence of the wild-type. The dots give the results of individual-based simulations (see §5.3). The vaccine-escape mutant is assumed to have the following phenotype (slow mutant in figures 6 and 2): $\alpha_m = 0.02, \beta_m = 7, \gamma_m = 2, \epsilon_m = 1, \mathcal{R}_m = 3.46$. Other parameter values: $v = \delta = 3 \times 10^{-4}, n = 10^6, \omega_v = \omega_r = 0.05, \alpha_w = 0.02, \beta_w = 10, \gamma_w = 2, \epsilon_w = 0.05, \mathcal{R}_w = 4.95$. The light grey area on the right-hand side indicates the speed above which the wild-type pathogen is expected to be driven to extinction ($v > v_c \approx 0.264$, see equation (3.2)).</figcaption></figure>

arise more frequently in infected hosts who are partially immunized [7–9]. Our model can be used to explore the consequences of this within-host evolution in vaccinated hosts (e.g. by taking $\theta_v > \theta_u$). A larger value of $\theta_v$ increases the overall rate of mutation (figure 3) but this effect is modulated by the fraction of the host population that is vaccinated. Consequently, when $\theta_v > \theta_u$, the speed of vaccination rollout can have a non-monotonic effect on the probability that a vaccine-escape mutation is introduced (see electronic supplementary material, figure SI.4). Indeed, when the rate of vaccination remains low, the enhancing effect of vaccination on the rate of introduction of new mutations can counteract the delaying effect of faster vaccination rollout discussed above. But the probability that a vaccine-escape mutation is introduced drops to very low levels when the rate of vaccination approaches the critical vaccination rate $v_c$.

The second step of adaptation starts when the vaccine-escape mutant has been introduced in the pathogen population. Will this new variant go extinct rapidly or will it start to invade? The answer to this question depends on the time at which the mutant is introduced. If the mutant is introduced when the population is not at an endemic equilibrium, the fate of the mutant depends on a time-varying infection rate which is driven by the fluctuations in the density of susceptible hosts. In our model, early introductions are likely to result in rapid extinction because there are simply not enough vaccinated hosts to favour the mutant over the wild-type. Moreover, we found that earlier introductions are likely to favour slower life-history strategies which are less prone to early extinction. If the introduction takes place later, when the system has reached a new endemic equilibrium, the fate of the
in the pathogen population. Will this new variant go extinct rapidly or will it start to invade? The answer to this question depends on the time at which the mutant is introduced. If the mutant is introduced when the population is not at an endemic equilibrium, the fate of the mutant depends on a time-varying birth rate which is driven by the fluctuations of the density of susceptible hosts. Earlier introductions are likely to result in rapid extinction because there are simply not enough vaccinated hosts to favour the mutant over the wild type. Moreover, we found that earlier introductions are likely to favour slower life-history strategies which are less prone to early extinction. If the introduction takes place later, when the system has reached a new endemic equilibrium, the fate of the mutant is solely governed by the effective per-generation ratio $R_m^*$ and does not depend on the life-history traits of the mutant. Slow and fast variants have equal probability to invade if they have the same $R_m^*$. Altogether, our results suggest that earlier arrival may not always facilitate invasion since the probability of invasion is limited by the time-varying epidemiological state of the host population.

The third step of adaptation starts as soon as the hosts infected by the vaccine-escape mutant are abundant and the effect of demographic stochasticity on the dynamics of this mutation becomes negligible. Our analysis attempts to better characterize the dynamics of the mutant after invasion using a combination of deterministic and analytical approximations. In principle, conditional on invasion, we can use the deterministic model (1) to describe the joint dynamics of the mutant and the wild type. In particular, the speed at which the vaccine-escape mutant spreads in the pathogen population can be well approximated by the deterministic model. This may be particularly useful to address the impact of various vaccination strategies on the speed of the spread of a vaccine-adapted variant [15]. In the present work we show that life-history traits of the vaccine-escape mutant drive the speed of its spread. Indeed, as pointed out before, the deterministic transient dynamics depends on the per-capita growth rate of the mutant $r_m$, not its per-generation reproduction ratio $R_m$ [13]. Transient dynamics may favour a fast and aggressive variant because this life-history strategy may be more competitive away from the endemic equilibrium. Yet, this explosive strategy may be risky for the pathogen if it leads to epidemiological fluctuations that result in a massive drop in the number of infections. The consequences of such fluctuations on the extinction risk of the variant
mutant is solely governed by the effective per-generation ratio $\mathcal{R}_m^*$ and does not depend on the life-history traits of the mutant. Slow and fast variants have equal probability to invade if they have the same $\mathcal{R}_m^*$. Altogether, our results suggest that earlier arrival may not always facilitate invasion since the probability of invasion is limited by the time-varying epidemiological state of the host population.

The third step of adaptation starts when the hosts infected by the vaccine-escape mutant are abundant and the effect of demographic stochasticity on the dynamics of this mutation becomes negligible. Our analysis attempts to better characterize the dynamics of the mutant after invasion using a combination of deterministic and stochastic approximations. In principle, conditional on invasion, we can use the deterministic model (2.4) to describe the joint dynamics of the mutant and the wild-type. In particular, the speed at which the vaccine-escape mutant spreads in the pathogen population can be well approximated by the deterministic model. This may be particularly useful to address the impact of various vaccination strategies on the speed of the spread of a vaccine-escape variant [12]. In the present work, we show that life-history traits of the vaccine-escape mutant drive the speed of its spread. Indeed, as pointed out before, the deterministic transient dynamics depends on the *per capita* growth rate of the mutant $r_m$, not its per-generation reproduction ratio $\mathcal{R}_m$ [45]. Transient dynamics may favour a fast and aggressive variant (i.e. faster increase in frequency of this variant) because this life-history strategy may be more competitive away from the endemic equilibrium. Yet, this explosive strategy may be risky for the pathogen if it leads to epidemiological fluctuations that result in a large drop in the number of infections. The consequences of such fluctuations on the extinction risk of the variant can be accounted for by a generalized birth-and-death process where the *per capita* growth rate of the mutant varies with time. Epidemiological fluctuations lead to a degradation of the future environment (i.e. depletion of the density of susceptible hosts) which results in an increased risk of extinction [46,50]. This effect has recently been analysed in a purely epidemiological model without vaccination [31]. In this simpler scenario, it is also relevant to make a distinction between early extinction (i.e. a *fizzle* in [31, fig. 1c]) and extinction after a successful invasion (i.e. an *epidemic burnout* in [31, fig. 1d]) and it is possible to use a similar hybrid semi-deterministic approach to obtain accurate analytical approximations for both events.

A comprehensive understanding of pathogen dynamics after vaccination relies on the use of a combination of theoretical tools to capture the interplay between stochastic and deterministic forces. Here, we use a hybrid numerical-analytical approach to account for the three successive steps that may eventually lead to the fixation of a vaccine-escape mutant. This theoretical framework is particularly suitable to explore the influence of different vaccination strategies on the risk of pathogen adaptation. In particular, we show that this risk drops to very low levels even when the speed of vaccination rollout is below the threshold value that may eventually lead to eradication (i.e. $v < v_c$). In other words, faster vaccination rollout makes sense even when eradication is infeasible, because faster rollout decreases both the number of cases and the likelihood of pathogen evolution. This conclusion is akin to the general prediction that the rate of pathogen adaptation should be maximized for intermediate immune pressure or for medium doses of chemotherapy at the within-host level [52,59-64]. Most of these earlier studies focused on evolutionary rescue scenarios where the wild-type is expected to be rapidly driven to extinction by human intervention. Our versatile theoretical framework, however, allows us to deal with a broader range of situations where the intervention is not expected to eradicate the wild-type pathogen. Accounting for the dynamics of the wild-type affects both the flux of mutation and the fate of these mutations. Note how our decomposition of the factors acting on the probability of adaptation (figure 6) provides a validation of the verbal argument often used in earlier studies to explain the higher rate of pathogen adaptation for intermediate levels of vaccination coverage of drug concentration [7,59-61].

The framework we have developed can be readily extended to explore many other situations. For instance, our model can be modified to explore the influence of temporal variations in the environment that could be driven by seasonality or by non-pharmaceutical interventions (NPIs). We explored a situation where the transmission rate of all variants is periodically reduced by a quantity $1 - c(t)$, where $c(t)$ is a measure of the intensity of NPIs. These periodic interventions affect both the flux of mutations and the probability that these mutations invade. In particular, NPIs lower the probability of mutant introduction through the reduction in the density of hosts infected by the wild-type (electronic supplementary material, figure SI.4). As a consequence, the probability of adaptation is reduced when vaccination is combined with periodic control measures. Hence, our approach helps to understand the interaction between vaccination and NPI discussed in earlier studies [14,65].

We have made several simplifying assumptions that need to be relaxed to confidently apply our findings to a broader range of pathogens (see electronic supplementary material, §6). First, one should study situations where the pathogen population has not reached an endemic equilibrium when vaccination starts to be applied. We carried out additional simulations showing that starting the vaccination rollout sooner (i.e. just after the start of the epidemic) tends to promote the probability of invasion of the escape mutant (electronic supplementary material, figure SI.5). Indeed, at the onset of the epidemic the density of susceptible hosts is higher (i.e. the birth rate of the infection is high relative to the endemic equilibrium) and the risk of early extinction of the mutant is reduced. Second, it is important to relax the assumption that natural immunity is perfect. We carried out additional simulations showing that when naturally immune hosts, like vaccinated hosts, can be reinfected the probability of invasion of the escape mutant increases (electronic supplementary material, figure SI.6). This effect is particularly strong just after the start of vaccination. Indeed, if naturally immune hosts are equivalent to vaccinated hosts, selection to escape immunity is present even before the start of vaccination and one may thus expect the speed of adaptation to be much faster. The vaccination strategy can, however, affect the rate of adaptation. In particular, we find that faster rates of vaccination always reduce the rate of adaptation via the reduction of the influx of escape mutants (electronic supplementary material, figure SI.7). Third, our model assumes that, once the infection takes place in a host, the vaccination status of the host does not affect the life cycle of the pathogen. Yet, if the vaccine reduces the transmission, the virulence and/or the duration of the infection we need to distinguish between infections of naive or vaccinated hosts. The effect of this heterogeneity falls beyond the scope of the present study. However, we note that this heterogeneity will require the use of multi-type branching processes to describe the invasion step of a rare variant [50] and it will require class-structured models to describe the deterministic dynamics of abundant variants [12]. This heterogeneity may
can be accounted for by a generalized birth-death process where the per-capita growth rate of the mutant varies with time. Epidemiological fluctuations lead to a degradation of the future environment (i.e., depletion of the density of susceptible hosts) which results in an increased risk of extinction [24, 6].

A comprehensive understanding of pathogen dynamics after vaccination relies on the use of a combination of theoretical tools to capture the interplay between stochastic and deterministic forces. Here, we use a hybrid numerical-analytical approach to account for the three successive steps that may eventually lead to the fixation of a vaccine-escape mutant. This theoretical framework is particularly suitable to explore the influence of different vaccination strategies on the risk of pathogen adaptation. In particular, we show that this risk drops to very low levels even when the speed of vaccination rollout is below the threshold value that may eventually lead to eradication (i.e., $\nu < \nu_c$). In other words, faster vaccination rollout makes sense even when eradication is infeasible, because faster rollout decreases both the number of cases and the likelihood of pathogen evolution. This conclusion is akin to the general prediction that the rate of pathogen adaptation should be maximized for intermediate immune pressure or for medium doses of chemotherapy at the within-host level [20, 40, 19, 1, 21, 11, 2]. Most of these earlier studies focused on evolutionary rescue scenarios where the wild type is expected to be rapidly driven to extinction by human intervention. Our versatile theoretical framework, however, allows us to deal with a broad range of situations where the intervention is not expected to eradicate the pathogen. Accounting for the dynamics of the wild type affects both the flux of mutation and the fate of these mutations.

The framework we have developed can be readily extended to explore may other situations. For instance, our model can be modified to explore the influence of temporal variations in the environment that could be driven by seasonality or by non-pharmaceutical interventions (NPIs). We explored a situation where the transmission rate of all variants is periodically reduced by a quantity $1-c(t)$, where $c(t)$ is a measure of the intensity of NPIs. These periodic interventions affect both the flux of mutations and the probability that these mutations invade (Figures S3 and S4). Notably, NPIs lower the probability of mutant introduction through the reduction in the density of hosts infected by the wild type. As a consequence, the probability of adaptation is reduced when
promote the evolutionary branching of the pathogen population and the long-term coexistence of different strains will require new theoretical developments to predict the fate of a new variant [66,67]. Another important extension of our model would be to study the effect of a diversity of vaccines in the host population. We did not explore this effect in the present study but this diversity of immune profiles among vaccinated hosts could slow down pathogen adaptation if escape from different vaccines requires distinct mutations [68–70].

Finally, it is important to recall that we focus here on a simplified scenario where we analyse the evolutionary epidemiology of an isolated population. In real-life situations the arrival time may depend more on the immigration of new variants from abroad than on local vaccination policies. The influence of migration remains to be investigated in spatially structured models where vaccination may vary among populations [71].

# 5. Methods

In this section, we present how extinction, invasion and fixation probabilities may be obtained under strong-selection assumptions when a mutant strain appears in a host-pathogen system that is away from its endemic equilibrium. Our essential tools are the deterministic ordinary differential equations (§2.1.1) and birth-and-death process approximations (§2.1.2). The former allows us to consider the situation when all strains are abundant, the latter when at least one strain is rare. We will limit ourselves to an informal treatment, presenting heuristic arguments and deferring rigorous proofs and sharp error bounds to a future treatment.

## 5.1. Approximating the fixation probability

Suppose that the mutant strain introduced at time $T_{int} = t_{int}$ successfully invades; we next consider the probability, $P_{fix}(t_{int})$, that the mutant will outcompete the wild-type and go to fixation. Fixation of the mutant occurs if it is still present when the wild-type strain disappears. If we let $T_{ext}^m$ and $T_{ext}^w$ be the extinction times of mutant and wild-type strains, respectively, then the probability of mutant fixation is $\mathbb{P}\{T_{ext}^w < T_{ext}^m\}$, which we may decompose as

$$ 
\begin{aligned}
\int_{t_{int}}^{\infty} \mathbb{P}\{T_{ext}^m > t\} \mathbb{P}\{T_{ext}^w \in [t, t + dt)\} &= - \int_{t_{int}}^{\infty} \mathbb{P}\{T_{ext}^m > t\} \frac{d}{dt} \mathbb{P}\{T_{ext}^w > t\} dt \\
&= - \int_{t_{int}}^{\infty} \mathbb{P}\{I_m^m(t) > 0\} \frac{d}{dt} (1 - \mathbb{P}\{I_w^w(t) = 0\}) dt \quad (5.1) \\
&= \int_{t_{int}}^{\infty} \mathbb{P}\{I_m^m(t) > 0\} \frac{d}{dt} \mathbb{P}\{I_w^w(t) = 0\} dt.
\end{aligned}
$$

We obtain estimates of $\mathbb{P}\{I_w^m(t) > 0\}$ and $\mathbb{P}\{I_m^m(t) > 0\}$ ($t > t_{int}$) by now approximating both mutant and wild-type strains by birth-and-death-processes $\tilde{I}_m(t)$ and $\tilde{I}_w(t)$ (see §2.1.2). The birth (i.e. infection) rates for the two types, $i = w, m$, are given by

$$ b_i(t) = \frac{\beta_i(X_u(t) + \epsilon_i X_v(t))}{N(t)} \quad (5.2) $$

and the death (i.e. death or recovery) rates are

$$ d_i = \delta + \alpha_i + \gamma_i \quad (5.3) $$

for $i \in \{w, m\}$.

As previously, we approximate the frequency of unvaccinated and vaccinated hosts by their deterministic counterparts:

$$ \frac{X_u^n(t)}{N^n(t)} \approx \frac{X_u(t)}{N(t)} \quad \text{and} \quad \frac{X_v^n(t)}{N^n(t)} \approx \frac{X_v(t)}{N(t)} $$

and compute the latter using the ordinary differential equation (2.4). Unlike previously, when we assumed that the mutant was rare and took $Y_m(t) \equiv 0$, we are now allowing the possibility that the mutant is abundant, and cannot neglect the effect of the mutant strain on $X_v$ and $X_u$. In particular, we need to take care in choosing the initial conditions of (2.4) to account for the fact that we consider the time of appearance of the first mutant that successfully invades. We are conditioning on the non-extinction of the mutant strain, and must account for the inherent variability in the time required to invade. This results in a random initial condition for the deterministic dynamics (see electronic supplementary material, §4, for details). In practice, we find that the randomness has negligible effect, but we must still take the conditioning into account. To do so, we first use (2.4) with $Y_m(0) = 0$ (so $Y_m(t) \equiv 0$ for $t > 0$) and initial conditions (electronic supplementary material, SI.1) to compute the epidemiological dynamics of the wild-type from time 0 up until the introduction of the mutant at time $t_{int}$. Then, at time $t_{int}$, we restart (2.4) with new initial conditions: we use the values $X_u(t_{int})$, $X_v(t_{int})$, $N(t_{int})$ and $Y_w(t_{int})$ computed assuming $Y_m(0) = 0$, and take

$$ Y_m(t_{int}) = \frac{1}{P_{inv}(t_{int})n} \quad (5.4) $$

(see electronic supplementary material, §4, for details). Crucially, the initial density of the mutant depends on the probability of successful invasion, $P_{inv}(t_{int})$, obtained above (3.7).
vaccination is combined with periodic control measures (Figure 6). Hence, our approach helps to understand the interaction between vaccination and NPI discussed in earlier studies [41, 31]. An important limitation of this study is that, for conceptual clarity, we have made several simplifying assumptions that would need to be relaxed to confidently apply our findings to the current SARS-CoV-2 pandemic. First, one should study situations where the pathogen population has not reached an endemic equilibrium when selection (e.g., vaccination or chemotherapy) starts to be applied. Second, it would be necessary to relax other assumptions made here, such as frequency-dependent transmission or life-long immunity after recovery. Accounting for natural recovery would require yet another class of imperfectly immune hosts. Similarly, another important extension of our model would be to study the effect of a diversity of vaccines in the host population. This diversity of immune profiles could slow down pathogen adaptation if the escape of different vaccines requires distinct mutations [45, 7, 35]. Finally, it is important to recall that we focus here on a simplified scenario where we analyse the evolutionary epidemiology of an isolated population. In real-life situations the arrival time may depend more on the immigration of new variants from abroad than on local vaccination policies. The influence of migration remains to be investigated in spatially structured models where vaccination may vary among populations [18].

# 5 Methods

In this section, we present how extinction, invasion and fixation probabilities may be obtained under strong-selection assumptions when a mutant strain appears in a host-pathogen system that is away from its endemic equilibrium. Our essential tools are the deterministic ordinary differential equations (1) and birth-and-death approximations, which we discuss below. The former allows us to consider the situation when all strains are abundant, the latter when at least one strain is rare. We will limit ourselves to an informal treatment, presenting heuristic arguments and deferring rigorous proofs and sharp error bounds to a future treatment. In the following we present a simple, yet versatile, hybrid (i.e., semi-deterministic) framework which allows us to approximate the probabilities associated with different steps of adaptation (mutation, invasion, fixation) via adding auxiliary equations describing stochastic phenomena, to the deterministic ordinary differential equations
Provided we use (2.4) with the appropriate initial conditions, then, as previously, the birth rates of both the wild-type and mutant strains are approximately deterministic, and we have

$$ \mathbb{P}\{I_i^n(t) > 0\} \approx \mathbb{P}\{\tilde{I}_i(t) > 0\}. \quad (5.5) $$

Under the branching assumption, the lines of descent of distinct infected individuals are independent, hence the probability that strain i vanishes by time t is the product of the probabilities that each line of descent vanishes:

$$ \mathbb{P}\{\tilde{I}_i(t) > 0\} = 1 - (1 - U_i(t|t_{\text{int}}))^{I_i^n(t_{\text{int}})} \approx 1 - (1 - U_i(t|t_{\text{int}}))^{nY_i(t_{\text{int}})}, \quad (5.6) $$

where, using [46],

$$ U_i(t|t_{\text{int}}) = \frac{1}{1 + \int_{t_{\text{int}}}^t d_i e^{-\int_{t_{\text{int}}}^s b_i(u) - d_i du} ds} \quad (5.7) $$

is the probability that an individual infected with strain i ∈ {w, m} present at time t_{int} has descendants alive at time t > t_{int}. In (5.6), we approximate the initial number of individuals infected with strain i using the frequencies obtained using (2.4) and (5.4):

$$ I_i^n(t_{\text{int}}) = nY_i^n(t_{\text{int}}) \approx nY_i(t_{\text{int}}). \quad (5.8) $$

Below in §5.2.3, we give a fast numerical method for computing $U_i(t|t_{\text{int}})$.

## 5.2. Auxiliary functions

In the following, we present a simple, yet versatile, hybrid (i.e. semi-deterministic and semi-numerical) framework which allows us to approximate the probabilities associated with different steps of adaptation (mutation, invasion, fixation). We do this by adding auxiliary equations describing stochastic phenomena to the deterministic ordinary differential equations describing the global population dynamics.

### 5.2.1. Introduction of the variant by mutation (step 1)

Recall $F_{\text{int}}(t) = 1 - e^{-\int_0^t \lambda_{\text{int}}(s) ds}$, (3.6), where $\lambda_{\text{int}}(t) = \theta_u Y_{uw}(t) + \theta_v Y_{vw}(t)$, (3.5). Rather than computing the integral—which would require that we compute $\lambda_{\text{int}}(s)$ (and thus $Y_{uw}(s)$ and $Y_{vw}(s)$) for every $s < t$— we observe that the cumulative hazard $\Lambda_{\text{int}}(t) = \int_0^t \lambda_{\text{int}}(s) ds$ can be computed by combining (2.4) with initial conditions (electronic supplementary material, SI.1) and the auxiliary differential equation

$$ \dot{\Lambda}_{\text{int}} = \lambda_{\text{int}} \quad (5.9) $$

with initial condition $\Lambda_{\text{int}}(0) = 0$. The use of this auxiliary equation reduces the computational effort by obtaining $\Lambda_{\text{int}}(t)$ simultaneously with $Y_{uw}(t)$ and $Y_{vw}(t)$ (as opposed to computing the latter two and then integrating).

### 5.2.2. Invasion of the variant (step 2)

In practice, the probability of mutant invasion (3.7) involves integrals that cannot be explicitly computed, and we must compute it numerically. To do so, we make use of one of the steps involved in computing $P_{\text{inv}}(t_{\text{int}})$ in [46]. There, it is shown that

$$ P_{\text{inv}}(t_{\text{int}}) = U_m(\infty|t_{\text{int}}) = \lim_{t \to \infty} U_m(t|t_{\text{int}}), $$

where

$$ U_m(t|t_{\text{int}}) = \mathbb{P}\{\tilde{I}_m(t) > 0 | \tilde{I}_m(t_{\text{int}}) = 1\} \quad (5.10) $$

is obtained via a pair of auxiliary functions

$$ \dot{U}_m = -d_m U_m V_m, \quad (5.11a) $$

$$ \dot{V}_m = (d_m - b_m(t))V_m - d_m V_m^2, \quad (5.11b) $$

with initial conditions

$$ U_m(t_{\text{int}}|t_{\text{int}}) = V_m(t_{\text{int}}|t_{\text{int}}) = 1 $$

(note that $\dot{U}_m$ and $\dot{V}_m$ denote the derivatives with respect to t). We compute $b_m(t)$, which depends on $X_u(t)$, $X_v(t)$ and $N(t)$ (see (2.8)), via (2.4). In practice, we cannot compute $U_m(\infty|t_{\text{int}})$; to obtain an approximation we approximate it by $U_m(t|t_{\text{int}})$ for the first value of t sufficiently large that

$$ |U_m(t + \Delta t|t_{\text{int}}) - U_m(t|t_{\text{int}})| $$

is less than our desired threshold of error, where $\Delta t$ is the step size in our numerical scheme.
374 describing the global population dynamics.

## 5.1 Before the introduction of a variant

376 We assume that vaccination starts after the monomorphic population of the wild-type pathogen
377 has reached its endemic equilibrium,

$$
\begin{align*}
S_U^0 &= \frac{\lambda(\delta + \gamma_w)}{\delta(\beta_w - \alpha_w)} \\
S_V^0 &= 0 \\
I_{Uw}^0 &= \frac{\lambda(R_w - 1)}{\beta_w - \alpha_w} \tag{11} \\
I_{Vw}^0 &= 0 \\
N^0 &= \frac{\lambda R_w (\delta + \gamma_w)}{\delta(\beta_w - \alpha_w)}
\end{align*}
$$

378 We then use the following ordinary differential equations to track the deterministic dynamics
379 of the wild-type pathogen using the endemic equilibrium before vaccination (11) as the initial
380 condition:

$$
\begin{align}
\dot{S}_U &= \lambda(1-p) + \omega S_V - \left( \beta_w \frac{I_{Uw} + I_{Vw}}{N} + \delta + \nu \right) S_U \nonumber \\
\dot{S}_V &= \lambda p + \nu S_U - \left( e_w \beta_w \frac{I_{Uw} + I_{Vw}}{N} + \delta + \omega \right) S_V \nonumber \\
\dot{I}_{Uw} &= \beta_w (I_{Uw} + I_{Vw}) \frac{S_U}{N} - (\delta + \alpha_w + \gamma_w) I_{Uw} \tag{12} \\
\dot{I}_{Vw} &= e_w \beta_w (I_{Uw} + I_{Vw}) \frac{S_V}{N} - (\delta + \alpha_w + \gamma_w) I_{Vw} \nonumber \\
\dot{N} &= \lambda - \delta N - \alpha_w (I_{Uw} + I_{Vw}). \nonumber
\end{align}
$$

Letting $I_w = I_{Uw} + I_{Vw}$, we get from (12):

$$
\dot{I}_w = \beta_w (S_U + e_w S_V) \frac{I_w}{N} - (\delta + \alpha_w + \gamma_w) I_w.
$$

A new endemic equilibrium will thus be approached after vaccination if and only if the growth rate
$r_w = \beta_w(S_U + e_w S_V)/N - (\delta + \alpha_w + \gamma_w)$ is positive, or equivalently if the effective per-generation
### 5.2.3. Fixation of the variant (step 3)

Here, we need two pairs of auxiliary equations to track the probabilities that some descendant of a wild-type or mutant individual that was present at $t_{int}$ is still alive at time $t$. Let

$$U_i(t|t_{int}) = \mathbb{P}\{\tilde{I}_i(t) > 0 | \tilde{I}_i = 1\}, \quad (5.12)$$

for $i \in \{\text{w}, \text{m}\}$. Exactly as in (5.11) above, these $U_i(t|t_{int})$ satisfy

$$\dot{U}_i = -d_i U_i V_i, \quad (5.13a)$$

$$\dot{V}_i = (d_i - b_i(t))V_i - d_i V_i^2, \quad (5.13b)$$

with

$$U_i(t_{int}|t_{int}) = V_i(t_{int}|t_{int}) = 1,$$

for $i \in \{\text{u}, \text{v}\}$.

To compute the probability of fixation, we first consider the probability that fixation occurs prior to time $t$, which is derived in exactly the same manner as (5.1):

$$\mathbb{P}\{\text{fixation prior to } t\} = \int_{t_{\text{int}}}^{t} \mathbb{P}\{I_{\text{m}}^{n}(s) > 0\} \frac{\mathrm{d}s}{\mathrm{d}s} \mathbb{P}\{I_{\text{w}}^{n}(s) = 0\} \mathrm{d}s.$$

Proceeding as in §5.1, we approximate the probabilities $\mathbb{P}\{I_{\text{m}}^{n}(s) > 0\}$ and $\mathbb{P}\{I_{\text{w}}^{n}(s) = 0\}$ by $\mathbb{P}\{\tilde{I}_{\text{m}}(s) > 0\}$ and $\mathbb{P}\{\tilde{I}_{\text{w}}(s) = 0\}$, and the initial number of hosts infected with the wild-type using the deterministic density, $I_{\text{w}}^{n}(t_{\text{int}}) \approx nY_{\text{w}}(t_{\text{int}})$. Then, using the branching property [24], the probability of fixation is approximately

$$U_{\text{fix}}(t|t_{\text{int}}) = \int_{t_{\text{int}}}^{t} \left(1 - (1 - U_{\text{m}}(s|t_{\text{int}}))^{nY_{\text{w}}(t_{\text{int}})}\right) \times \left(nY_{\text{w}}(t_{\text{int}})(-\dot{U}_{\text{w}}(s|t_{\text{int}}))(1 - U_{\text{w}}(s|t_{\text{int}}))^{nY_{\text{w}}(t_{\text{int}})-1}\right) \mathrm{d}s. \quad (5.14)$$

Differentiating yields the following auxiliary equation for $U_{\text{fix}}(t, t_{\text{int}})$:

$$\dot{U}_{\text{fix}} = nY_{\text{w}}(t_{\text{int}})(\delta + \alpha_{\text{w}} + \gamma_{\text{w}})U_{\text{w}}V_{\text{w}}(1 - U_{\text{w}})^{nY_{\text{w}}(t_{\text{int}})-1}U_{\text{m}}, \quad (5.15)$$

with initial condition $U_{\text{fix}}(t_{\text{int}}|t_{\text{int}}) = 0$. We estimate the fixation probability as

$$P_{\text{fix}}(t_{\text{inv}}) = \lim_{t \to \infty} U_{\text{fix}}(t|t_{\text{int}}), \quad (5.16)$$

approximating the limit at infinity as we did for $P_{\text{inv}}(t_{\text{int}})$ in §5.2.2 above.

### 5.2.4. The overall risk of pathogen adaptation

We numerically compute the cumulative density function $F_{\text{inv}}(t) = \mathbb{P}\{T_{\text{inv}} \le t\}$ of the first arrival time $T_{\text{inv}}$ of a vaccine-escape mutant that successfully invades (3.10) analogously to $F_{\text{int}}$ (§5.2.1), using the auxiliary equation

$$\Lambda_{\text{inv}} = \lambda_{\text{int}} P_{\text{inv}} \quad (5.17)$$

with initial condition $\Lambda_{\text{inv}}(0) = 0$. We compute $Y_{\text{uw}}(t)$ and $Y_{\text{vw}}(t)$—and thus $\lambda_{\text{int}}(t)$—using (2.4) with initial conditions (electronic supplementary material, SI.1).

### 5.3. Stochastic simulations

We carried out stochastic simulations to check the validity of our results. We developed an individual-based simulation program for the Markov process described in table 1, using the parameter values given in table 2. In order to match the assumptions used in our analysis, we started the simulation when the system was at its endemic equilibrium before vaccination. We then introduced a single host infected with the mutant pathogen at a time $t_{\text{int}}$ after the start of vaccination and let the simulation run until one of the pathogen variants (the wild-type or the mutant) went extinct. If the wild-type goes extinct first we recorded this run as a 'mutant fixation event'. We ran 1000 replicates for each set of parameters and plotted the proportion of runs that led to mutant fixation in figure 5. We also used our simulations to confirm our prediction on the speed of viral adaptation in figure 6. In this scenario, we allowed the vaccine-escape variant to be introduced by mutation from the wild-type genotype. We carried out 1000 simulations and monitored (i) the frequency of the escape mutant at different points in time after the start of vaccination (figure 6a) and (ii) the number of introduction events by mutation (figure 6b). We also used this simulation approach to go beyond the scenarios used in our analysis to check the robustness of some of our results in electronic supplementary material, §6.

**Ethics.** This work did not require ethical approval from a human subject or animal welfare committee.

**Data accessibility.** The simulation code used to carry out stochastic simulations has been deposited on Zenodo [72].
reproduction ratio

$$R_w^e = R_w \left( \frac{S_U}{N} + e_w \frac{S_V}{N} \right)$$

is larger than 1, when $S_U$ and $S_V$ are taking their stationary values $\tilde{S}_U$ and $\tilde{S}_V$ in the absence of infection. Computing these values (see SI, Section 1) shows that $R_w^e > 1$ if and only if $e_w R_w > 1$ or $\nu > \nu_c$, where

$$\nu_c = \frac{R_w(\delta(1 - (1 - e_w)p) + \omega) - (\delta + \omega)}{1 - R_w e_w} \quad (13)$$

Thus we see that if $e_w R_w < 1$ and the speed of the vaccination rollout is higher than the critical value $\nu_c$ the wild type will be driven to extinction deterministically.

For values of the vaccination rollout $\nu$ smaller than this threshold $\nu_c$, or when $e_w R_w > 1$, the wild-type may also become extinct because of demographic stochasticity. We can neglect this possibility because the timescale of stochastic extinction from abundances of the order of $n$ is much larger than those of the processes under consideration.

## 5.2 Introduction of the variant by mutation (step 1)

As indicated above, we use a time-inhomogeneous Poisson point process to model the influx of new mutations. The per capita rate of mutation is assumed to be constant through time but whether or not a mutant will escape extinction within a host may depend on the type of host. Indeed, a vaccine-escape mutation may have a higher probability to escape within-host extinction in vaccinated hosts. We account for this effect by making a distinction between $\theta_U$ and $\theta_V$. If vaccine-escape mutations are more likely to escape extinction in vaccinated hosts we expect $\theta_V > \theta_U$. In other words, $\theta_V/\theta_U - 1$ is a measure of the within-host fitness advantage of the vaccine-escape mutant in vaccinated hosts (they are assumed to have the same within-host fitness in naïve hosts).

We can compute the probability that some of the vaccine-escape mutations are present as standing variation before the start of vaccination. When the resident population has reached its endemic equilibrium ($S_U^0, 0, I_{Uw}^0, 0, N^0$), the number of mutants is approximated by a birth-and-death process with immigration, with birth rate $b_m^0 = \beta_m \frac{S_U^0}{N^0}$ and death rate $d_m = \delta + \alpha_m + \gamma_m$, and the “immigration” is actually mutations arising in the resident population, which occur at rate
Electronic supplementary material is available online [73].

**Declaration of AI use.** We have not used AI-assisted technologies in creating this article.

**Authors' contributions.** S.G.: conceptualization, formal analysis, funding acquisition, investigation, supervision, visualization, writing—original draft, writing—review and editing; A.L.: conceptualization, formal analysis, writing—review and editing; V.M.: software, visualization, writing—review and editing; T.D.: conceptualization, writing—review and editing; T.L.P.: conceptualization, formal analysis, investigation, supervision, writing—original draft, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** S.G. acknowledges support from the CNRS PEPS 2022 grant 'VaxDurable' and the ANR VIRADAPT (ANR-24-CE35-5671).

**Acknowledgements.** We thank reviewers for constructive comments.

## References

1. Gandon S, Mackinnon M, Nee S, Read A. 2003 Imperfect vaccination: some epidemiological and evolutionary consequences. *Proc. R. Soc. Lond. B* **270**, 1129–1136. (doi:10.1098/rspb.2003.2370)

2. McLean AR. 1995 After the honeymoon in Measles Control. *Lancet* **345**, 272. (doi:10.1016/s0140-6736(95)90272-4)

3. Gandon S, Day T. 2008 Evidences of parasite evolution after vaccination. *Vaccine* **26**, C4–C7. (doi:10.1016/j.vaccine.2008.02.007)

4. Kennedy DA, Read AF. 2017 Why does drug resistance readily evolve but vaccine resistance does not? *Proc. R. Soc. Lond. B* **284**, 20162562. (doi:10.1098/rspb.2016.2562)

5. Kennedy DA, Read AF. 2018 Why the evolution of vaccine resistance is less of a concern than the evolution of drug resistance. *Proc. Natl Acad. Sci. USA* **115**, 12878–12886. (doi:10.1073/pnas.1717159115)

6. Restif O, Grenfell BT. 2007 Vaccination and the dynamics of immune evasion. *J. R. Soc. Interface* **4**, 143–153. (doi:10.1098/rsif.2006.0167)

7. Saad-Roy CM et al. 2021 Epidemiological and evolutionary considerations of SARS-CoV-2 vaccine dosing regimes. *Science* **372**, 363–370. (doi:10.1126/science.abg8663)

8. Cobey S, Larremore DB, Grad YH, Lipsitch M. 2021 Concerns about SARS-CoV-2 evolution should not hold back efforts to expand vaccination. *Nat. Rev. Immunol.* **21**, 330–335. (doi:10.1038/s41577-021-00544-9)

9. Day T, Kennedy DA, Read AF, Gandon S. 2022 Pathogen evolution during vaccination campaigns. *PLoS Biol.* **20**, e3001804. (doi:10.1371/journal.pbio.3001804)

10. Rouzine IM, Rozhnova G. 2023 Evolutionary implications of SARS-CoV-2 vaccination for the future design of vaccination strategies. *Commun. Med.* **3**, 86. (doi:10.1038/s43856-023-00320-x)

11. Gandon S. 2018 Evolutionary epidemiology theory of vaccination. In *Pertussis: epidemiology, immunology, and evolution* (eds P Rohani, SV Scarpino), p. 113. Oxford, UK: Oxford University Press.

12. Gandon S, Lion S. 2022 Targeted vaccination and the speed of SARS-CoV-2 adaptation. *Proc. Natl Acad. Sci. USA* **119**, e2110666119. (doi:10.1073/pnas.2110666119)

13. Thompson RN, Hill EM, Gog JR. 2021 SARS-CoV-2 incidence and vaccine escape. *Lancet Infect. Dis.* **21**, 913–914. (doi:10.1016/S1473-3099(21)00202-4)

14. Rella SA, Kulikova YA, Dermitzakis ET, Kondrashov FA. 2021 Rates of SARS-CoV-2 transmission and vaccination impact the fate of vaccine-resistant strains. *Sci. Rep.* **11**, 10. (doi:10.1038/s41598-021-95025-3)

15. Kogan O, Khasin M, Meerson B, Schneider D, Myers CR. 2014 Two-strain competition in quasineutral stochastic disease dynamics. *Phys. Rev. E* **90**, 042149. (doi:10.1103/PhysRevE.90.042149)

16. Humplik J, Hill AL, Nowak MA. 2014 Evolutionary dynamics of infectious diseases in finite populations. *J. Theor. Biol.* **360**, 149–162. (doi:10.1016/j.jtbi.2014.06.039)

17. Parsons TL, Lambert A, Day T, Gandon S. 2018 Pathogen evolution in finite populations: slow and steady spreads the best. *J. R. Soc. Interface* **15**, 20180135. (doi:10.1098/rsif.2018.0135)

18. Sanjuán R, Nebot MR, Chirico N, Mansky LM, Belshaw R. 2010 Viral mutation rates. *J. Virol.* **84**, 9733–9748. (doi:10.1128/JVI.00694-10)

19. Day T, Proulx SR. 2004 A general theory for the evolutionary dynamics of virulence. *Am. Nat.* **163**, E40–E63. (doi:10.1086/382548)

20. Day T, Gandon S. 2007 Applying population-genetic models in theoretical evolutionary epidemiology. *Ecol. Lett.* **10**, 876–888. (doi:10.1111/j.1461-0248.2007.01091.x)

21. Hethcote HW. 1989 Three basic epidemiological models. In *Applied mathematical ecology* (eds SA Levin, TG Hallam, LJ Gross), pp. 119–144. Berlin, Germany: Springer. (doi:10.1007/978-3-642-61317-3_5)

22. Kermack W, McKendrick A. 1927 A contribution to the mathematical theory of epidemics. *Proc. R. Soc. Lond. A* **115**, 700–721. (doi:10.1098/rspa.1927.0118)

23. Gillespie D. 1977 Exact stochastic simulation of coupled chemical reactions. *J. Phys. Chem* **81**, 2340–2361. (doi:10.1021/j100540a008)

24. Bianchi F, Mascipinto S, Stefanizzi P, De Nitto S, Germinario C, Tafuri S. 2020 Long-term immunogenicity after measles vaccine vs. wild infection: an Italian retrospective cohort study. *Hum. Vaccines Immunother.* **17**, 2078–2084. (doi:10.1080/21645515.2020.1871296)

25. Kim JH, Liepkalns J, Reber AJ, Lu X, Music N, Jacob J, Sambhara S. 2016 Prior infection with influenza virus but not vaccination leaves a long-term immunological imprint that intensifies the protective efficacy of antigenically drifted vaccine strains. *Vaccine* **34**, 495–502. (doi:10.1016/j.vaccine.2015.11.077)

26. Chen YQ et al. 2018 Influenza infection in humans induces broadly cross-reactive and protective neuraminidase-reactive antibodies. *Cell* **173**, 417–429. (doi:10.1016/j.cell.2018.03.030)

27. Viboud C et al. 2020 Beyond clinical trials: evolutionary and epidemiological considerations for development of a universal influenza vaccine. *PLoS Pathog.* **16**, e1008583. (doi:10.1371/journal.ppat.1008583)

28. Greaney AJ, Loes AN, Gentles LE, Crawford KHD, Starr TN, Malone KD, Chu HY, Bloom JD. 2021 Antibodies elicited by mRNA-1273 vaccination bind more broadly to the receptor binding domain than do those from SARS-CoV-2 infection. *Sci. Transl. Med.* **13**, i9915. (doi:10.1126/scitranslmed.abi9915)

29. McCallum H. 2001 How should pathogen transmission be modelled? *Trends Ecol. Evol.* **16**, 295–300. (doi:10.1016/s0169-5347(01)02144-9)

30. Parsons TL. 2012 Asymptotic analysis of some stochastic models from population dynamics and population genetics. PhD thesis, [Canada]: University of Toronto.

31. Parsons TL, Bolker BM, Dushoff J, Earn DJD. 2024 The probability of epidemic burnout in the stochastic SIR model with vital dynamics. *Proc. Natl Acad. Sci. USA* **121**, e2313708120. (doi:10.1073/pnas.2313708120)

32. Bartlett MS. 1956 Deterministic and stochastic models for recurrent epidemics (ed. J Neyman). In *Proc. 3rd Berkely Symp. on Mathematical Statistics and Probability*, pp. 81–109. Berkeley, CA: University of California Press. (doi:10.1525/9780520350717-007)

33. Ethier SN, Kurtz TG. *Markov processes: characterization and convergence*, pp. 803–817. Hoboken, NJ: Wiley Interscience.
404 $\mu_m = \theta_U I_{Uw}^0$. Because we assume that in a fully naïve host population vaccine-escape mutations 405 carry a fitness cost relative to the resident strain, we have $b_m^0 < d_m$. The number of mutants 406 thus approximately follows a subcritical birth-and-death process with immigration, which is known 407 to converge in distribution as $t \to \infty$ to a negative binomial stationary distribution [30]. The 408 probability that there are $k$ infected individuals hosting the mutant pathogen at time $t = 0$ is thus:

$$p_k = \binom{k+r-1}{k} (1-R_m^0)^r (R_m^0)^k \quad (14)$$

409 where $r = \mu_m/b_m^0$ and $R_m^0 = b_m^0/d_m$. Hence the the expected number of vaccine-escape mutants 410 already present at the start of vaccination is

$$\frac{\mu_m}{d_m - b_m} \quad (15)$$

411 This result is analogous to the classical result that the expected frequency of deleterious mutations 412 is of the form $\mu/s$ where $\mu$ is the rate of mutation towards deleterious mutants and $s$ is the fitness 413 cost of those deleterious mutants.

414 We can also compute the probability that no mutant is present at the start of vaccination:

$$p_0 = (1 - R_m^0)^r \quad (16)$$

415 When either $R_m^0$ or $r$ is small, $p_0 \approx 1$, and we can neglect the presence of preexisting mutants. 416 Otherwise, we need to account for the possibility that one or more mutants are present at time 417 $t = 0$, which we discuss in Section 5.5 below.

418 We now assume that there is zero mutant present at the start of vaccination. We are interested 419 in the law of the first time $t_m$ at which a mutant appears. Because $\theta_U I_{Uw}(t) + \theta_V I_{Vw}(t)$ is the flux 420 of vaccine-escape mutants from the wild-type population, by the exponential formula for Poisson 421 point processes, we have [27]:

$$\mathbb{P}\{t_m > t\} = e^{-\int_0^t (\theta_U I_{Uw}(s) + \theta_V I_{Vw}(s))ds} \quad (17)$$
34. Andersson H, Britton T. 2000 *Stochastic epidemic models and their statistical analysis*. New York, NY: Springer.

35. Jagers P. 1992 Stabilities and instabilities in population dynamics. *J. Appl. Probab.* **29**, 770–780. (doi:10.1017/S0021900200043667)

36. Kribs-Zaleta CM, Velasco-Hernández JX. 2000 A simple vaccination model with multiple endemic states. *Math. Biosci.* **164**, 183–201. (doi:10.1016/S0025-5564(00)00003-1)

37. Hadeler KP, van den Driessche P. 1997 Backward bifurcation in epidemic control. *Math. Biosci.* **146**, 15–35. (doi:10.1016/S0025-5564(97)00027-8)

38. Allen L. 2010 *An introduction to stochastic processes with applications to biology*. Boca Raton, FL: CR press.

39. Hamza K, Jagers P, Klebaner FC. 2016 On the establishment, persistence, and inevitable extinction of populations. *J. Math. Biol.* **72**, 797–820. (doi:10.1007/s00285-015-0903-2)

40. van Herwaarden OA, Grasman J. 1995 Stochastic epidemics: major outbreaks and the duration of the endemic period. *J. Math. Biol.* **33**, 581–601. (doi:10.1007/BF00298644)

41. Andersson H, Djehiche B. 1998 A threshold limit theorem for the stochastic logistic epidemic. *J. Appl. Probab.* **35**, 662–670.

42. Parsons TL. 2018 Invasion probabilities, hitting times, and some fluctuation theory for the stochastic logistic process. *J. Math. Biol.* **77**, 1193–1231. (doi:10.1007/s00285-018-1250-x)

43. Ewens WJ. 1979 *Mathematical population genetics*. Berlin, Germany: Springer-Verlag.

44. Cox D, Isham V. 1980 *Point processes*. vol. 12. Boca Raton, FL: CRC Press.

45. Gandon S, Day T. 2007 The evolutionary epidemiology of vaccination. *J. R. Soc. Interface* **4**, 803–817. (doi:10.1098/rsif.2006.0207)

46. Kendall DG. 1948 On the generalized 'birth-and-death' process. *Ann. Math. Statist.* **19**, 1–15. (doi:10.1214/aoms/1177730285)

47. Geritz SAH. 2005 Resident-invader dynamics and the coexistence of similar strategies. *J. Math. Biol.* **50**, 67–82. (doi:10.1007/s00285-004-0280-8)

48. Cai Y, Geritz SAH. 2020 Resident-invader dynamics of similar strategies in fluctuating environments. *J. Math. Biol.* **81**, 907–959. (doi:10.1007/s00285-020-01532-8)

49. Priklopil T, Lehmann L. 2020 Invasion implies substitution in ecological communities with class-structured populations. *Theor. Popul. Biol.* **134**, 36–52. (doi:10.1016/j.tpb.2020.04.004)

50. Carmona P, Gandon S. 2020 Winter is coming: pathogen emergence in seasonal environments. *PLoS Comput. Biol.* **16**, e1007954. (doi:10.1371/journal.pcbi.1007954)

51. Martin G, Aguilée R, Ramsayer J, Kaltz O, Ronce O. 2013 The probability of evolutionary rescue: towards a quantitative comparison between theory and evolution experiments. *Phil. Trans. R. Soc. B* **368**, 20120088. (doi:10.1098/rstb.2012.0088)

52. Alexander HK, Martin G, Martin OY, Bonhoeffer S. 2014 Evolutionary rescue: linking theory for conservation and medicine. *Evol. Appl.* **7**, 1161–1179. (doi:10.1111/eva.12221)

53. Bell G. 2017 Evolutionary rescue. *Annu. Rev. Ecol. Evol. Syst.* **48**, 605–627. (doi:10.1146/annurev-ecolsys-110316-023011)

54. Rausch JW, Capoferri AA, Katusiime MG, Patro SC, Kearney MF. 2020 Low genetic diversity may be an Achilles heel of SARS-CoV-2. *Proc. Natl Acad. Sci. USA* **117**, 24614–24616. (doi:10.1073/pnas.2017726117)

55. Davies NG et al. 2021 Estimated transmissibility and impact of SARS-CoV-2 lineage B.1.1.7 in England. *Science* **372**, eabg3055. (doi:10.1126/science.abg3055)

56. Volz E et al. 2021 Assessing transmissibility of SARS-CoV-2 lineage B.1.1.7 in England. *Nature* **593**, 266–269. (doi:10.1038/s41586-021-03470-x)

57. Wang P et al. 2021 Antibody resistance of SARS-CoV-2 variants B.1.351 and B.1.1.7. *Nature* **593**, 130–135. (doi:10.1038/s41586-021-03398-2)

58. Paton RS, Overton CE, Ward T. 2022 The rapid replacement of the SARS-CoV-2 Delta variant by Omicron (B.1.1.529) in England. *Sci. Transl. Med.* **14**, eabo5395. (doi:10.1126/scitranslmed.abo5395)

59. Grenfell BT, Pybus OG, Gog JR, Wood JLN, Daly JM, Mumford JA, Holmes EC. 2004 Unifying the epidemiological and evolutionary dynamics of pathogens. *Science* **303**, 327–332. (doi:10.1126/science.1090727)

60. Read AF, Day T, Huijben S. 2011 The evolution of drug resistance and the curious orthodoxy of aggressive chemotherapy. *Proc. Natl Acad. Sci. USA* **108**, 10871–10877. (doi:10.1073/pnas.1100299108)

61. Kouyos RD et al. 2014 The path of least resistance: aggressive or moderate treatment? *Proc. R. Soc. Lond. B* **281**, 20140566. (doi:10.1098/rspb.2014.0566)

62. Alexander HK, Bonhoeffer S. 2012 Pre-existence and emergence of drug resistance in a generalized model of intra-host viral dynamics. *Epidemics* **4**, 187–202. (doi:10.1016/j.epidem.2012.10.001)

63. Hartfield M, Alizon S. 2015 Within-host stochastic emergence dynamics of immune-escape mutants. *PLoS Comput. Biol.* **11**, e1004149. (doi:10.1371/journal.pcbi.1004149)

64. Day T, Read AF. 2016 Does high-dose antimicrobial chemotherapy prevent the evolution of resistance? *PLoS Comput. Biol.* **12**, e1004689. (doi:10.1371/journal.pcbi.1004689)

65. Lobinska G, Pauzner A, Traulsen A, Pilpel Y, Nowak MA. 2022 Evolution of resistance to COVID-19 vaccination with dynamic social distancing. *Nat. Hum. Behav.* **6**, 193–206. (doi:10.1038/s41562-021-01281-8)

66. Park SW, Cobey S, Metcalf CJE, Levine JM, Grenfell BT. 2024 Predicting pathogen mutual invasibility and co-circulation. *Science* **386**, 175–179. (doi:10.1126/science.adq0072)

67. Lion S, Sasaki A, Boots M. 2023 Extending eco-evolutionary theory with oligomorphic dynamics. *Ecol. Lett.* **26**, S22–S46. (doi:10.1111/ele.14183)

68. van Houte S et al. 2016 The diversity-generating benefits of a prokaryotic adaptive immune system. *Nature* **532**, 385–388. (doi:10.1038/nature17436)

69. Chabas H, Lion S, Nicot A, Meaden S, van Houte S, Moineau S, Wahl LM, Westra ER, Gandon S. 2018 Evolutionary emergence of infectious diseases in heterogeneous host populations. *PLoS Biol.* **16**, e2006738. (doi:10.1371/journal.pbio.2006738)

70. McLeod D, Wahl L, Mideo N. 2021 Mosaic vaccination: how distributing different vaccines across a population could improve epidemic control. *Evol. Lett.* **5**, 458–471. (doi:10.1002/evl3.252)

71. Gerrish PJ et al. 2021 How unequal vaccine distribution promotes the evolution of vaccine escape. SSRN. (doi:10.2139/ssrn.3827009)

72. Gandon S, Lambert A, Voinson M, Day T, Parsons TL. 2025 The speed of vaccination rollout and the risk of pathogen adaptation. Zenodo. (doi:10.5281/zenodo.12655541)

73. Gandon S, Lambert A, Voinson M, Day T, Parsons TL. 2025 Supplementary material from: The speed of vaccination rollout and the risk of pathogen adaptation. Figshare. (doi:10.6084/m9.figshare.c.7828935)
We can numerically compute the probability density function $f_m$ of the first arrival time $t_m$ of a vaccine-escape mutant using

$$f_m(t) = \dot{F}_m(t)e^{-F_m(t)}, \qquad (18)$$

where $F_m(t)$ is given by the auxiliary equation

$$\dot{F}_m = \theta_U I_{Uw}(t) + \theta_V I_{Vw}(t) \qquad (19)$$

with initial condition $F_m(0) = 0$, while we compute $I_{Uw}(t)$ and $I_{Vw}(t)$ using (12). The use of this auxiliary equation reduces computational effort by obtaining $F_m(t)$ simultaneously with $I_{Uw}(t)$ and $I_{Vw}(t)$ (as opposed to computing the latter and then integrating).

## 5.3 Invasion of the variant (step 2)

Suppose that a mutant strain appears at time $t_m \ge 0$ in a single infected host, that is, with density $I_m(0) = \frac{1}{n}$, (which is effectively zero as $n$ becomes large). Then, (2) yields $I_m(t_m) = 0$ for all $t$, whereas the dynamics of the system follows (12). This differential equation approximation does not mean that the mutant is absent, but simply not sufficient in numbers to be visible at the coarse resolution and short time scale upon which (12) is applicable.

Then we combine (12) with a birth-and-death process approximation, $\tilde{I}_m(t)$, to the number of individuals infected with the mutant strain at time $t$ after the mutant arrival time $t_m$, $I_m^n(t)$. We approximate the rate of new infections, $\frac{\beta_m(S_U^n(t)+e_mS_V^n(t))}{N^n(t)}$ by replacing the stochastic quantities $S_U^n(t)$, $S_V^n(t)$ and $N^n(t)$ by their deterministic approximations, giving the time-dependent birth rate

$$b_m(t) = \frac{\beta_m(S_U(t) + e_mS_V(t))}{N(t)}$$

where $S_U(t)$, $S_V(t)$ and $N(t)$ are determined via the deterministic system (12). Each death in the birth-and-death process corresponds to the removal of a susceptible, which occurs by host death or recovery at combined rate $d_m = \delta + \alpha_m + \gamma_m$. See §8.2 in the Supplementary Information of [36] for a rigorous justification.

The so-called “merciless dichotomy” [23] tells us that the time-inhomogeneous birth-and-death
process started with one individual at time $t_m$ either goes extinct, or grows without bound (i.e. invades) with probability (see [24] and SI, Section 2)

$$P_m^{t_m} = \frac{1}{1 + \int_{t_m}^{\infty} d_m e^{-\int_{t_m}^{t} b_m(s) - d_m \, ds}} dt$$

Thus, either the mutant strain vanishes, or the number infected with the mutant strain will eventually be of the order of $n$ individuals, after which we can use (1) to compute the densities of both wild-type and mutant strains.

In practice, we can compute the probability of mutant invasion when the mutant is introduced at time $t_m$ using $P_m^{t_m} = U_m^{t_m}(\infty)$ where $U_m^{t_m}(\infty)$ is obtained via the pair of auxiliary functions $U_m^{t_m}$ and $V_m^{t_m}$ [24] defined as follows: $U_m^{t_m}(t) = \mathbb{P}\{\tilde{I}_m(t) \neq 0\}$ and $V_m^{t_m}(t) = \mathbb{P}\{\tilde{I}_m(t) = 1 | \tilde{I}_m(t) \neq 0\}$, We then have

$$\dot{U}_m^{t_m} = -d_m U_m^{t_m} V_m^{t_m} \qquad (20)$$

$$\dot{V}_m^{t_m} = (d_m - b_m(t))V_m^{t_m} - d_m(V_m^{t_m})^2, \qquad (21)$$

where $U_m^{t_m}(t_m) = V_m^{t_m}(t_m) = 1$ and we compute $b_m(t)$, $S_U(t)$, $S_V(t)$ and $N(t)$, via (12). In practice, we cannot compute $U_m^{t_m}(\infty)$; to obtain an approximation we evaluate $U_m^{t_m}(t)$ for sufficiently large $t$ that $|U_m^{t_m}(t + \Delta t) - U_m^{t_m}(t)|$ is less than our desired threshold of error.

Note that several variants can arise and fail to invade before finally a lucky variant manages to invade. We can use the probability of invasion $P_m^t$ of a variant introduced at time $t$ to characterize the distribution of the first time $t_i$ at which a mutant is introduced that successfully invades. By the thinning property of Poisson point processes, we have [27]:

$$\mathbb{P}\{t_i > t\} = e^{-\int_0^t (\theta_U I_{Uw}(s) + \theta_V I_{Vw}(s)) P_m^s ds} \qquad (22)$$

where $\theta_U I_{Uw}(t) + \theta_V I_{Vw}(t)$ is the flux of vaccine-escape mutants from the wild-type population. We compute numerically the probability density function $g_m$ of the first arrival time $t_i$ of a vaccine-
458 escape mutant that successfully invades using

$$g_m(t) = \dot{G}_m(t)e^{-G_m(t)}, \qquad (23)$$

459 where $G_m(t)$ is given by the auxiliary equation

$$\dot{G}_m(t) = (\theta_U I_{Uw}(t) + \theta_V I_{Vw}(t)) P_m^t \qquad (24)$$

460 with initial condition $G_m(0) = 0$ and computing $I_{Uw}(t)$ and $I_{Vw}(t)$ using (12).

Compare (22) with (17) and note that the probability that no vaccine-escape mutant will ever arise is

$$\mathbb{P}\{t_m = \infty\} = e^{-\int_0^\infty (\theta_U I_{Uw}(t) + \theta_V I_{Vw}(t)) dt}.$$

In contrast, the probability that no vaccine-escape mutant will ever invade is the larger probability

$$\mathbb{P}\{t_i = \infty\} = e^{-\int_0^\infty (\theta_U I_{Uw}(t) + \theta_V I_{Vw}(t)) P_m^t dt}.$$

Note that $P_m^t$ converges as $t \to \infty$ to $P_m^* = 1 - 1/R_m^*$ which is nonzero, so that

$$\mathbb{P}\{t_m = \infty\} = 0 \Leftrightarrow \int_0^\infty (\theta_U I_{Uw}(t) + \theta_V I_{Vw}(t)) dt = \infty \Leftrightarrow \mathbb{P}\{t_i = \infty\} = 0,$$

461 that is, the probability of adaptation is 1 if and only if $t \mapsto (\theta_U I_{Uw}(t) + \theta_V I_{Vw}(t))$ is not integrable.

462 In other words, the probability of adaptation is 1 in the limit $t \to \infty$ when the wild type is not driven to extinction by vaccination (i.e. $\nu < \nu_c$) which implies that there is an uninterrupted flux of 463 mutation producing vaccine-escape variants. One of these mutants will eventually escape extinction 464 and invade. Yet, the time needed for a successful variant to appear may be very long and we focus 465 in the main text on $A_m^t$ the probability of adaptation before time $t$ (equation (10) and Figure 6). 466
### 5.4 Fixation of the variant (step 3)

Suppose now that the mutant strain successfully invades; we next consider the probability $P_{fix}$ that the mutant will outcompete the wild type and go to fixation. Fixation of the mutant occurs if it is still present when the wild-type strain disappears. If we let $T_m$ and $T_w$ be the extinction times of mutant and wild-type strains, the probability of mutant fixation is thus $\mathbb{P}\{T_w < T_m\}$ which we may decompose as

$$ 
\begin{aligned}
\int_{t_m}^{\infty} \mathbb{P}\{T_m > t\}\mathbb{P}\{T_w \in dt\} &= -\int_{t_m}^{\infty} \mathbb{P}\{T_m > t\} \frac{d}{dt}\mathbb{P}\{T_w > t\} dt \\
&= -\int_{t_m}^{\infty} \mathbb{P}\{I_m^n(t) > 0\} \frac{d}{dt}(1 - \mathbb{P}\{I_w^n(t) = 0\}) dt \quad (25) \\
&= \int_{t_m}^{\infty} \mathbb{P}\{I_m^n(t) > 0\} \frac{d}{dt}\mathbb{P}\{I_w^n(t) = 0\} dt.
\end{aligned}
$$

We again obtain estimates of $\mathbb{P}\{I_w^n(t) > 0\}$ and $\mathbb{P}\{I_m^n(t) > 0\}$ using the fact that conditional on $S_U(t)$, $S_V(t)$ and $N(t)$, $(\tilde{I}_w(t), \tilde{I}_m(t))$ follows a time-inhomogeneous, two-type birth-and-death process, where the birth rates for the two types, $i = w, m$, are given by

$$ b_i(t) = \frac{\beta_i(S_U(t) + e_i S_V(t))}{N(t)} $$

and the death rates are $d_i = \delta + \alpha_i + \gamma_i$. The birth rates vary with time due to the epidemiological perturbations following the start of vaccination and in particular, to the feedback of mutant invasion on $S_V$ and $S_U$. To quantify these epidemiological perturbations, we now approximate the density of susceptibles and total host density by the values of $S_U(t)$, $S_V(t)$ and $N(t)$ obtained from the full deterministic system (1), which accounts for the presence of the mutant by replacing $I_m$ with its *expected value*.

We need to take care in choosing the initial conditions of (1) to account for the fact that we consider the time of appearance of the first mutant that successfully invades and so are conditioning on the non-extinction of the mutant strain, and for the inherent variability in the time required to invade; this results in a random initial condition for the deterministic dynamics (see Supplementary Information §4 for details). In practice, we find that the randomness has negligible effect, but
we must still take the conditioning into account. To do so, we first use (12) to compute the epidemiological dynamics of the wild type before the introduction of the mutant at time $t_m$. Then, at time $t_m$, we use (1) where $S_U(t_m)$, $S_V(t_m)$, $N(t_m)$ and $I_w(t_m)$ are obtained from (12) and take $I_m(t_m) = \frac{1}{(1-P_m^{t_m})n}$ (see Supplementary Information §4). Crucially, the initial density of the mutant depends on the probability of successful invasion of the mutant $P_m^{t_m}$ obtained above.

Provided we use (1) with the appropriate initial conditions as previously, the birth rates of both the wild-type and mutant strains are approximately deterministic, and from [24], we have:

$$ \mathbb{P}\{I_i^n(t) > 0\} \approx \mathbb{P}\{\tilde{I}_i(t) > 0\} = 1 - (1 - U_i^{t_m}(t))^{I_i^n(t_m)}, \quad (26) $$

where

$$ U_i^{t_m}(t) = \frac{1}{1 + \int_{t_m}^{t} d_i e^{-\int_{t_m}^{u} b_i(s) - d_i ds} du} \quad (27) $$

is the probability that an individual infected with strain $i$ ($i = m, w$) present at $t_m$ has descendants alive at time $t$. Under the branching assumption, the lines of descent of distinct infected individuals are independent, hence the probability that strain $i$ vanishes by time $t$ is the product of the probabilities that each line of descent vanishes, $(1 - U_i^{t_m}(t))^{I_i^n(t_m)}$.

In practice, we need two pairs of auxiliary equations to track the probability of extinction of both the wild type and the mutant

$$ \dot{U}_w^{t_m} = -d_w U_w^{t_m} V_w^{t_m} \quad (28) $$

$$ \dot{V}_w^{t_m} = (d_w - b_w(t))V_w^{t_m} - d_w(V_w^{t_m})^2 \quad (29) $$

$$ \dot{U}_m^{t_m} = -d_m U_m^{t_m} V_m^{t_m} \quad (30) $$

$$ \dot{V}_m^{t_m} = (d_m - b_m(t))V_m^{t_m} - d_m(V_m^{t_m})^2, \quad (31) $$

with $U_w^{t_m}(t_m) = V_w^{t_m}(t_m) = U_m^{t_m}(t_m) = V_m^{t_m}(t_m) = 1$. To compute the probability of fixation, we first consider the probability of fixation prior to time $t$, which is derived in exactly the same manner as (25).

$$ U^{t_m}(t) = \int_{t_m}^{t} \mathbb{P}\{I_m^n(s) > 0\} \frac{d}{ds} \mathbb{P}\{I_w^n(s) = 0\} ds, $$
502 Using (26) to approximate the probabilities $P\{I_m^n(s) > 0\}$ and $P\{I_w^n(s) = 0\}$, we get

$$U^{t_m}(t) = \int_{t_m}^{t} [1 - (1 - U_m^{t_m}(s))^{I_m^n(t_m)}][I_w^n(t_m)(-\dot{U}_w^{t_m})(1 - U_w^{t_m}(s))^{I_w^n(t_m)-1}]ds.$$

503 Differentiating this and taking $I_w^n(t_m) = nI_w(t_m)$ and $I_m^n(t_m) = 1$ yields the following auxiliary equation for $U^{t_m}(t)$:

$$\dot{U}^{t_m} = nI_w(t_m)(\delta + \alpha_w + \gamma_w)U_w^{t_m}V_w^{t_m}(1 - U_w^{t_m})^{nI_w(t_m)-1}U_m^{t_m}, \quad (32)$$

505 with initial condition $U^{t_m}(t_m) = 0$. We estimate the fixation probability as $P_{fix}^{t_m} = U^{t_m}(\infty)$ as 506 above.

## 5.5 Invasion and Fixation with Standing Variation

508 We now come back to the probability of adaptation from standing variation at time $t = 0$, using the probability $p_k$ that there are $k$ mutants present at time $t = 0$ and the estimates for the probabilities 509 of invasion and fixation $P_m^{t_m}$ and $P_{fix}^{t_m}$, of a single mutant arriving at time $t_m$, taking $t_m = 0$. Recall 510 that $p_k$ is negative binomial with success probability $R_m^0 = b_m^0/d_m$ and $r = \mu_m/b_m^0$ failures. Under 511 the branching process approximation, the chain of infections started by each mutation will go 512 extinct independently with probability $1 - P_m^0$. The probability of invasion is then the probability 513 that at least one line survives, $1 - (1 - P_m^0)^k$. Summing this over all possible values of $k$ gives us 514 the invasion probability from standing variation, 515

$$\sum_{k=1}^{\infty} p_k (1 - (1 - P_m^0)^k) = 1 - p_0 - \sum_{k=1}^{\infty} p_k (1 - P_m^0)^k = 1 - \sum_{k=0}^{\infty} p_k (1 - P_m^0)^k.$$

516 Recalling that the probability generating function for the number of mutants at time $t = 0$ is

$$\sum_{k=0}^{\infty} p_k z^k = \left( \frac{1 - R_m^0}{1 - R_m^0 z} \right)^r,$$

517 which converges provided $|z| < \frac{1}{R_m^0}$, we see that the probability of invasion from standard variation 518 is $1 - \left( \frac{1 - R_m^0}{1 - R_m^0(1 - P_m^0)} \right)^r$.
Similarly, if there are $k$ mutations at time $t = 0$, (26) gives us that $\mathbb{P}\{I_m^n(t) > 0 | I_m(0) = k\} \approx 1 - (1 - U_m^0(t))^k$, so proceeding as above, we find that the probability that the mutant is still present at time $t$, assuming that at least one individual was present at time $t = 0$ is approximately $1 - \left(\frac{1 - R_m^0}{1 - R_m^0(1 - U_m^0(t))}\right)^r$. Substituting this for $\mathbb{P}\{I_m^n(t) > 0\}$ in (25) and differentiating as above gives us an auxiliary equation analogous to (32) for $U^0(t)$, the probability that the mutant fixes starting from standing variation:

$$ 
\begin{aligned}
\dot{U}^0 &= \mathbb{P}\{I_m^n(s) > 0\} \frac{d}{dt} \mathbb{P}\{I_w^n(s) = 0\} \\
&= nI_w(0)(\delta + \alpha_w + \gamma_w)U_w^0 V_w^0 (1 - U_w^0)^{nI_w(0)-1} \left(1 - \left(\frac{1 - R_m^0}{1 - R_m^0(1 - U_m^0)}\right)^r\right).
\end{aligned}
$$

As previously, we obtain $P_{fix}^0 = \dot{U}^0(\infty)$ by choosing $t$ sufficiently large that $\dot{U}^0(t)$ equilibrates.

## 5.6 Stochastic simulations

We carried out stochastic simulations to check the validity of our results (see **Figure S1**). We used an individual-based simulation program for the Markov process described in **Table S1**. We start the simulation when the system is at its endemic equilibrium before vaccination given by equation (11) in section 5.1. Then we introduce a single host infected with the mutant pathogen at a time $t_m$ after the start of vaccination and we let the simulation run until one of the pathogen variants (the wild-type or the mutant) goes extinct. If the wild-type goes extinct first we record this run as a “mutant fixation event”. We ran 1000 replicates for each set of parameters and **Figure S1** plots the proportion of runs that led to mutant fixation. Simulation code is available upon request and will be deposited on zenodo after acceptance of the manuscript.

**Acknowledgements:** SG acknowledges support from the CNRS PEPS 2022 grant "VaxDurable".

## References

[1] Helen K Alexander and Sebastian Bonhoeffer. Pre-existence and emergence of drug resistance in a generalized model of intra-host viral dynamics. *Epidemics*, 4(4):187–202, 2012.
535

536

537

[2] Helen K Alexander, Guillaume Martin, Oliver Y Martin, and Sebastian Bonhoeffer. Evolutionary rescue: linking theory for conservation and medicine. *Evolutionary applications*, 7(10):1161-1179, 2014.

538

539

[3] H. Andersson and T. Britton. *Stochastic Epidemic Models and their Statistical Analysis*, volume 151 of *Lecture notes in statistics*. Springer, New York, 2000.

540

541

[4] Graham Bell. Evolutionary rescue. *Annual Review of Ecology, Evolution, and Systematics*, 48:605-627, 2017.

542

543

[5] Yuhua Cai and Stefan A. H. Geritz. Resident-invader dynamics of similar strategies in fluctuating environments. 2020.

544

545

[6] Philippe Carmona and Sylvain Gandon. Winter is coming: Pathogen emergence in seasonal environments. *PLoS computational biology*, 16(7):e1007954, 2020.

546

547

548

[7] Hélène Chabas, Sébastien Lion, Antoine Nicot, Sean Meaden, Stineke van Houte, Sylvain Moineau, Lindi M Wahl, Edze R Westra, and Sylvain Gandon. Evolutionary emergence of infectious diseases in heterogeneous host populations. *PLoS biology*, 16(9):e2006738, 2018.

549

550

551

[8] Sarah Cobey, Daniel B Larremore, Yonatan H Grad, and Marc Lipsitch. Concerns about sars-cov-2 evolution should not hold back efforts to expand vaccination. *Nature Reviews Immunology*, pages 1-6, 2021.

552

553

554

555

[9] Nicholas G Davies, Sam Abbott, Rosanna C Barnard, Christopher I Jarvis, Adam J Kucharski, James D Munday, Carl AB Pearson, Timothy W Russell, Damien C Tully, Alex D Washburne, et al. Estimated transmissibility and impact of sars-cov-2 lineage b. 1.1. 7 in england. *Science*, 372(6538), 2021.

556

557

[10] Troy Day, David A Kennedy, Andrew F Read, and Sylvain Gandon. The evolutionary epidemiology of pathogens during vaccination campaigns. *arXiv preprint arXiv:2109.13680*, 2021.

558

559

[11] Troy Day and Andrew F Read. Does high-dose antimicrobial chemotherapy prevent the evolution of resistance? *PLoS computational biology*, 12(1):e1004689, 2016.
560

561

[12] Saad-Roy et al. Epidemiological and evolutionary considerations of sars-cov-2 vaccine dosing regimes. medRxiv, 2020.

562

563

[13] S. Gandon and T. Day. The evolutionary epidemiology of vaccination. *J. Royal Soc. Interface*, 4(16):803–817, 2007.

564

565

[14] Sylvain Gandon and Troy Day. Evidences of parasite evolution after vaccination. *Vaccine*, 26:C4–C7, 2008.

566

567

[15] Sylvain Gandon and Sébastien Lion. Targeted vaccination and the speed of sars-cov-2 adaptation. medRxiv, 2021.

568

569

570

[16] Sylvain Gandon, Margaret Mackinnon, Sean Nee, and Andrew Read. Imperfect vaccination: some epidemiological and evolutionary consequences. *Proceedings of the Royal Society of London. Series B: Biological Sciences*, 270(1520):1129–1136, 2003.

571

572

[17] S. A. H. Geritz. Resident-invader dynamics and the coexistence of similar strategies. *J. Math. Biol.*, 50(1):67–82, 2005.

573

574

575

[18] Philip J Gerrish, Fernando Saldana, Benjamin Galeota-Sprung, Alexandre Colato, Erika E Rodriguez, and Jorge X Velasco-Hernandez. How unequal vaccine distribution promotes the evolution of vaccine escape. Available at SSRN 3827009, 2021.

576

577

578

[19] Bryan Grenfell, Nicole Mideo, Laura C Pollitt, Andrew F Read, David L Smith, Claire Standley, Nina Wale, Geoffrey Chi-Johnston, Ted Cohen, Troy Day, et al. The path of least resistance: aggressive or moderate. 2014.

579

580

581

[20] Bryan T Grenfell, Oliver G Pybus, Julia R Gog, James LN Wood, Janet M Daly, Jenny A Mumford, and Edward C Holmes. Unifying the epidemiological and evolutionary dynamics of pathogens. *science*, 303(5656):327–332, 2004.

582

583

[21] Matthew Hartfield and Samuel Alizon. Within-host stochastic emergence dynamics of immune-escape mutants. *PLoS computational biology*, 11(3):e1004149, 2015.
584

[22] Jan Humplik, Alison L Hill, and Martin A Nowak. Evolutionary dynamics of infectious diseases in finite populations. *Journal of theoretical biology*, 360:149-162, 2014.

585

586

[23] P. Jagers. Stabilities and instabilities in population dynamics. *J. Appl. Prob.*, pages 770-780, 1992.

587

588

[24] D. G. Kendall. On the generalized “Birth-and-Death” process. *Ann. Math. Stat.*, 19(1):1-15, 1948.

589

590

[25] David A Kennedy and Andrew F Read. Why does drug resistance readily evolve but vaccine resistance does not? *Proceedings of the Royal Society B: Biological Sciences*, 284(1851):20162562, 2017.

591

592

593

[26] David A Kennedy and Andrew F Read. Why the evolution of vaccine resistance is less of a concern than the evolution of drug resistance. *Proceedings of the National Academy of Sciences*, 115(51):12878-12886, 2018.

594

595

596

[27] John Frank Charles Kingman. *Poisson processes*, volume 3. Clarendon Press, 1992.

597

[28] Oleg Kogan, Michael Khasin, Baruch Meerson, David Schneider, and Christopher R Myers. Two-strain competition in quasineutral stochastic disease dynamics. *Physical Review E*, 90(4):042149, 2014.

598

599

600

[29] T. G. Kurtz. *Approximation of Population Processes*, volume 36. Society for Industrial and Applied Mathematics, Philadelphia, 1981.

601

602

[30] Amaury Lambert. Population dynamics and random genealogies. *Stochastic Models*, 24(sup1):45-163, 2008.

603

604

[31] Gabriela Lobinska, Ady Pauzner, Arne Traulsen, Yitzhak Pilpel, and Martin A Nowak. Evolution of resistance to covid-19 vaccination with dynamic social distancing. *Nature human behaviour*, 6(2):193-206, 2022.

605

606

607

[32] Guillaume Martin, Robin Aguilée, Johan Ramsayer, Oliver Kaltz, and Ophélie Ronce. The probability of evolutionary rescue: towards a quantitative comparison between theory and

608
evolution experiments. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 368(1610):20120088, 2013.

[33] Hamish McCallum, Nigel Barlow, and Jim Hone. How should pathogen transmission be modelled? *Trends in ecology & evolution*, 16(6):295-300, 2001.

[34] AngelaR McLean. After the honeymoon in measles control. *The Lancet*, 345(8945):272, 1995.

[35] David V McLeod, Lindi M Wahl, and Nicole Mideo. Mosaic vaccination: how distributing different vaccines across a population could improve epidemic control. *bioRxiv*, pages 2020-11, 2021.

[36] T. L. Parsons, A. Lambert, T. Day, and S. Gandon. Pathogen evolution in finite populations: slow and steady spreads the best. *J. Royal Soc. Interface*, 15(147), 2018.

[37] Robert S Paton, Christopher E Overton, and Thomas Ward. The rapid replacement of the delta variant by omicron (b. 1.1. 529) in england. *Sci. Transl. Med.*, page eabo5395, 2022.

[38] Tadeas Priklopil and Laurent Lehmann. Invasion implies substitution in ecological communities with class-structured populations. 134:36-52, 2020.

[39] Jason W Rausch, Adam A Capoferri, Mary Grace Katusiime, Sean C Patro, and Mary F Kearney. Low genetic diversity may be an achilles heel of sars-cov-2. *Proceedings of the National Academy of Sciences*, 117(40):24614-24616, 2020.

[40] Andrew F Read, Troy Day, and Silvie Huijben. The evolution of drug resistance and the curious orthodoxy of aggressive chemotherapy. *Proceedings of the National Academy of Sciences*, 108(Supplement 2):10871-10877, 2011.

[41] Simon A Rella, Yuliya A Kulikova, Emmanouil T Dermitzakis, and Fyodor A Kondrashov. Rates of sars-cov-2 transmission and vaccination impact the fate of vaccine-resistant strains. *Scientific Reports*, 11(1):1-10, 2021.

[42] Olivier Restif and Bryan T Grenfell. Vaccination and the dynamics of immune evasion. *Journal of the Royal Society Interface*, 4(12):143-153, 2007.
634

635

[43] Rafael Sanjuán, Miguel R Nebot, Nicola Chirico, Louis M Mansky, and Robert Belshaw. Viral mutation rates. *Journal of virology*, 84(19):9733-9748, 2010.

636

637

[44] Robin N Thompson, Edward M Hill, and Julia R Gog. Sars-cov-2 incidence and vaccine escape. *The Lancet Infectious Diseases*, 2021.

638

639

640

641

[45] Stineke van Houte, Alice KE Ekroth, Jenny M Broniewski, Hélène Chabas, Ben Ashby, Joseph Bondy-Denomy, Sylvain Gandon, Mike Boots, Steve Paterson, Angus Buckling, et al. The diversity-generating benefits of a prokaryotic adaptive immune system. *Nature*, 532(7599):385-388, 2016.

642

643

644

[46] Erik Volz, Swapnil Mishra, Meera Chand, Jeffrey C Barrett, Robert Johnson, Lily Geidelberg, Wes R Hinsley, Daniel J Laydon, Gavin Dabrera, Áine O'Toole, et al. Assessing transmissibility of sars-cov-2 lineage b. 1.1. 7 in england. *Nature*, 593(7858):266-269, 2021.

645

646

647

[47] Pengfei Wang, Manoj S Nair, Lihong Liu, Sho Iketani, Yang Luo, Yicheng Guo, Maple Wang, Jian Yu, Baoshan Zhang, Peter D Kwong, et al. Antibody resistance of sars-cov-2 variants b. 1.351 and b. 1.1. 7. *Nature*, pages 1-6, 2021.
<figure><img src="image_2.png" /><figcaption>**Figure 1: Graphical representation of the different evolutionary epidemiology outcomes after vaccination.** The density of the wild-type pathogen is indicated in light blue and the dynamics of the mutant in orange. Each panel describes the temporal dynamics of the epidemics after the start of vaccination: (a) eradication of the wild-type pathogen, (b) new endemic equilibrium of the wild-type population after damped oscillations (with no introduction of the vaccine-adapted mutant), (c) early extinction of the vaccine-adapted mutant after its introduction by mutation, (d) invasion of the vaccine-adapted mutant followed by its extinction, (e) invasion of the vaccine-adapted mutant and long-term coexistence with the wild type in a new endemic equilibrium after damped oscillations, (f) invasion and fixation of the vaccine-adapted mutant (extinction of the wild type). Note that (b), (c) and (d) result to the same endemic equilibrium (wild-type population) after damped oscillations. The vertical dashed line (black) indicates the start of vaccination. For simplicity we consider that vaccination starts after the wild-type population has reached an endemic equilibrium. The horizontal dashed line indicates the “stochastic threshold” above which one may consider that the deterministic model provides a very good approximation of the dynamics and we can neglect the effect of demographic stochasticity.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2:** A schematic representation of the model. Naive and uninfected hosts ($S_U$ hosts) are introduced at a rate $\lambda$ and are vaccinated with probability $p$ at birth and at rate $\nu$ later on. Immunization induced by the vaccine wanes at rate $\omega$. Uninfected hosts ($S_U$ and $S_V$) die at a rate $\delta$ while infected hosts ($I_{Ui}$ and $I_{Vi}$) die at a rate $d_i = \delta + \alpha_i$, where $i$ refers to the virus genotype: the wild-type ($i = w$) or the escape mutant ($i = m$). The rate of infection of naive hosts by the genotype $i$ is $h_i = \beta_i(I_{Ui} + I_{Vi})$, where $\beta_i$ is the transmission rate of the genotype $i$. Vaccination reduces the force of infection and $e_i$ refers to the ability of the genotype $i$ to escape the immunity triggered by vaccination (we assume $e_m > e_w$). A host infected by parasite genotype $i$ recovers from the infection at rate $\gamma_i$ and yields lifelong and perfect immunity ($R$ hosts). $N = S_U + S_V + \sum_i (I_{Ui} + I_{Vi}) + R$ is the total host population density.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Arrival time $t_m$ of the vaccine-escape mutant</figcaption></figure>

**Figure 3: Faster vaccine rollout delays the arrival time of the mutant.** We plot the probability density function $f_m(t)$ of the arrival time of the mutant for different speeds of vaccination rollout: $\nu = 0.1$ (top), 0.2 (middle) and 0.3 (bottom). We contrast a scenario where $\theta_V = \theta_U$ (dashed line), and $\theta_V = 10 \times \theta_U$ (full line). Other parameter values: $\theta_U = 1$, $\lambda = 0.01$, $\delta = 0.01$, $\gamma = 1$, $\beta_w = 20$, $\alpha_w = 1$, $e_w = 0.03$, $\omega = 0.05$, $p = 0.1$. For these parameter values the critical rate of vaccination $\nu_c$ above which the wild-type pathogen is driven to extinction is $\nu_c \approx 0.54$ (see equation (5)).
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Probability of invasion of the vaccine-escape mutant increases with $t_m$. We plot the probability invasion $P_m^{t_m}$ of a slow (green) and a fast (red) vaccine-escape mutant for different speeds of vaccination rollout: $v = 0.1$ (top), 0.2 (middle) and 0.3 (bottom). The slow mutant: $e_m = 1, \beta_m = 10, \alpha_m = 1$. The fast mutant: $e_m = 1, \beta_m = 30, \alpha_m = 5.02$. The probability of invasion $P_m^*$ in the limit $t_m \to \infty$ (see equation (9)) is indicated with the dashed black line. Other parameter values: $\lambda = 0.01, \delta = 0.01, \gamma = 1, \beta_w = 20, \alpha_w = 1, e_w = 0.03, \omega = 0.05, p = 0.1$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5:** Probability of fixation of the vaccine-escape mutant may be low when *t*<sub>m</sub> is large. We plot the probability of fixation of a slow (green) and a fast (red) vaccine-escape mutant for different speeds of vaccination rollout: *v* = 0.1 (top), 0.2 (middle) and 0.3 (bottom). The slow mutant: *e*<sub>m</sub> = 1, *β*<sub>m</sub> = 10, *α*<sub>m</sub> = 1. The fast mutant: *e*<sub>m</sub> = 1, *β*<sub>m</sub> = 30, *α*<sub>m</sub> = 5.02. For comparison with **Figure 4** we plot the probability of invasion *P*<sub>m</sub><sup>*t*<sub>m</sub></sup> with dashed colored lines and its asymptotic value *P*<sub>m</sub><sup>*</sup> with a dashed black line. Other parameter values: λ = 0.01, δ = 0.01, γ = 1, *β*<sub>w</sub> = 20, *α*<sub>w</sub> = 1, *e*<sub>w</sub> = 0.03, *ω* = 0.05, *p* = 0.1, *n* = 10<sup>8</sup>.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6: Probability of adaptation and the speed of vaccination rollout.** We plot the probability of adaptation $A_m^t$ of a vaccine-escape mutant against the speed of vaccination rollout for different values of time $t$ (black line). The dashed line indicates the probability of adaptation when we impose periodic fluctuations in $c(t)$ which measures the intensity of Non-Pharmaceutical Interventions (NPIs) that reduce the transmission rate of all pathogens by $1-c(t)$. Here we use a square wave function for $c(t)$ that fluctuates between 0.2 and 0 with a period $T = 200$. These NPIs affect both the flux of mutations (**Figure S3**) and the probability of invasion (**Figure S4**). The light gray area on the right-hand-side indicates the speed above which the wild-type pathogen is expected to be driven to extinction ($\nu > \nu_c$). In the absence of NPIs the critical rate of vaccination $\nu_c$ above which the wild-type pathogen is driven to extinction is $\nu_c \approx 0.77$ (see equation (5)). Parameter values: $\lambda = 0.01$, $\delta = 0.01$, $\gamma = 1$, $\beta_w = 20$, $\alpha_w = 1$, $e_w = 0.03$, $\beta_m = 10$, $\alpha_m = 1$, $e_m = 1$, $\omega = 0.05$, $p = 0$, $\theta_V = \theta_U = 0.1$.</figcaption></figure>