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.
1 Introduction
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 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.
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 the 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.
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 Sn, In and Rn, 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 Nn = Sn + In + Rn. 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 λ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 ν 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 and four classes of infected individuals
. We assume that the virulence αi (the mortality rate induced by the infection), the transmission βi (the production rate of new infections), and the recovery γi (the rate at which the host clears the infection) are fully governed by the pathogen genotype (i = w or m). A fourth trait ei ∈ [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. 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.
Naive and uninfected hosts (SU hosts) are introduced at a rate λ and are vaccinated with probability p at birth and at rate ν later on. Immunization induced by the vaccine wanes at rate ω. Uninfected hosts (SU and SV) die at a rate δ while infected hosts (IUi and IV i) die at a rate di = δ + α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 hi = βi(IUi + IV i), where βi is the transmission rate of the genotype i. Vaccination reduces the force of infection and ei refers to the ability of the genotype i to escape the immunity triggered by vaccination (we assume em > ew). A host infected by parasite genotype i recovers from the infection at rate γi and yields lifelong and perfect immunity (R hosts). N = SU + SV + Σ i(IUi + IV i) + R is the total host population density.
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 . 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:
It is also convenient to track the dynamics of the total density of hosts infected with the same strain i, Ii := IUi + IV i, which yields:
The ability of the strains to grow is given by the sign of the growth rate ri. Note that this growth rate depends on the four different traits of the pathogen: αi, βi, γi, ei. But this growth rate depends also on the densities SU and SV, which vary with t, the time since the start of vaccination (i.e., vaccination starts at t = 0). The coefficient of selection sm on the mutant strain relative to the wild type is:
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.
3.2 Pathogen eradication and vaccination threshold
The ability of the wild-type pathogen population to grow can be measured by its effective pergeneration reproduction ratio which is given by:
where
. 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
. Ultimately, vaccination can even lead to the eradication of the wild-type pathogen (Figure 1a) either when the vaccine is sufficiently efficient (ewRw > 1) or when the vaccination coverage is sufficiently high [34, 16]. The deterministic model (1) can be used to identify the threshold νc of the speed of vaccination rollout above which the wild-type pathogen can be driven to extinction (see Methods and SI Section 1):
As expected, better vaccines (i.e., lower values of ew and ω) yield lower threshold values for the speed of vaccination. Imperfect vaccines (i.e., higher values of ew and ω), 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 ν < ν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 ν < ν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 ν < ν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 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 and
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 θU /n if unvaccinated and θV /n if vaccinated. In addition, we assume that θU and θV are small enough that within-host clonal interference is unlikely, and that θV ≥ θ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
so that the probability density fm(t) of the arrival time t = tm of the vaccine-escape mutant is approximated by:
In other words, the time tm 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 fm for different values of rollout speed ν 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 θ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.
We plot the probability density function fm(t) of the arrival time of the mutant for different speeds of vaccination rollout: ν = 0.1 (top), 0.2 (middle) and 0.3 (bottom). We contrast a scenario where θV = θU (dashed line), and θV = 10 × θU (full line). Other parameter values: θU = 1, λ = 0.01, δ = 0.01, γ = 1, βw = 20, αw = 1, ew = 0.03, ω = 0.05, p = 0.1. For these parameter values the critical rate of vaccination νc above which the wild-type pathogen is driven to extinction is νc ≈ 0.54 (see equation (5)).
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 and constant death rate dm = δ + αm + γm, introduced at time tm, successfully invades (see [24] and SI, Section 2) is:
Plotting the probability of invasion vs. tm 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 SV (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)).
We plot the probability invasion of a slow (green) and a fast (red) vaccine-escape mutant for different speeds of vaccination rollout: ν = 0.1 (top), 0.2 (middle) and 0.3 (bottom). The slow mutant: em = 1, βm = 10, αm = 1. The fast mutant: em = 1, βm = 30, αm = 5.02. The probability of invasion
in the limit tm → ∞ (see equation (9)) is indicated with the dashed black line. Other parameter values: λ = 0.01, δ = 0.01, γ = 1, βw = 20, αw = 1, ew = 0.03, ω = 0.05, p = 0.1.
Taking tm → ∞ 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 and
, 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:
By definition, at the endemic equilibrium set by the wild-type pathogen we have
Hence, a necessary condition for the mutant to invade this equilibrium is
, 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
. This 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):
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
).
Interestingly, unlike , the probability
of mutant invasion at time tm given in (7) is not governed solely by Ri, but rather depends on the life-history traits of the mutants. For instance, assume that two vaccine-escape mutants have the same values of Rm and em 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 , 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 principle may be violated if the phenotype of the vaccine-escape mutant is very different than the phenotype of the wild type.
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 ri, not the per-generation reproduction ratio Ri [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 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 tm 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).
We plot the probability of fixation of a slow (green) and a fast (red) vaccine-escape mutant for different speeds of vaccination rollout: ν = 0.1 (top), 0.2 (middle) and 0.3 (bottom). The slow mutant: em = 1, βm = 10, αm = 1. The fast mutant: em = 1, βm = 30, αm = 5.02. For comparison with Figure 4 we plot the probability of invasion with dashed colored lines and its asymptotic value
with a dashed black line. Other parameter values: λ = 0.01, δ = 0.01, γ = 1, βw = 20, αw = 1, ew = 0.03, ω = 0.05, p = 0.1, n = 108.
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):
We plot the probability of adaptation 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 (ν > νc). In the absence of NPIs the critical rate of vaccination νc above which the wild-type pathogen is driven to extinction is νc ≈ 0.77 (see equation (5)). Parameter values: λ = 0.01, δ = 0.01, γ = 1, βw = 20, αw = 1, ew = 0.03, βm = 10, αm = 1, em = 1, ω = 0.05, p = 0, θV = θU = 0.1.
As expected, when ν < νc this probability of adaptation goes to 1 when t → ∞. Indeed, when vaccination cannot eradicate the wild type, a vaccine-adapted variant will eventually appear by mutation and invade. When ν > ν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 ν < ν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, ν < νc, that may ultimately lead to eradication. When the speed of vaccination is just below ν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 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 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 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
. 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 rm, not its per-generation reproduction ratio Rm [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 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., ν < ν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 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 describing the global population dynamics.
5.1 Before the introduction of a variant
We assume that vaccination starts after the monomorphic population of the wild-type pathogen has reached its endemic equilibrium,
We then use the following ordinary differential equations to track the deterministic dynamics of the wild-type pathogen using the endemic equilibrium before vaccination (11) as the initial condition:
Letting Iw = IUw + IV w, we get from (12):
A new endemic equilibrium will thus be approached after vaccination if and only if the growth rate rw = βw(SU + ewSV)/N − (δ + αw + γw) is positive, or equivalently if the effective per-generation reproduction ratio
is larger than 1, when SU and SV are taking their stationary values
and
in the absence of infection. Computing these values (see SI, Section 1) shows that
if and only if ewRw > 1 or ν > νc, where
Thus we see that if ewRw < 1 and the speed of the vaccination rollout is higher than the critical value νc the wild type will be driven to extinction deterministically.
For values of the vaccination rollout ν smaller than this threshold νc, or when ewRw > 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 θU and θV. If vaccine-escape mutations are more likely to escape extinction in vaccinated hosts we expect θV > θU. In other words, θV /θ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 , the number of mutants is approximated by a birth-and-death process with immigration, with birth rate
and death rate dm = δ + αm + γm, and the “immigration” is actually mutations arising in the resident population, which occur at rate
. Because we assume that in a fully naïve host population vaccine-escape mutations carry a fitness cost relative to the resident strain, we have
. The number of mutants thus approximately follows a subcritical birth-and-death process with immigration, which is known to converge in distribution as t → ∞ to a negative binomial stationary distribution [30]. The probability that there are k infected individuals hosting the mutant pathogen at time t = 0 is thus:
where
and
. Hence the the expected number of vaccine-escape mutants already present at the start of vaccination is
This result is analogous to the classical result that the expected frequency of deleterious mutations is of the form µ/s where µ is the rate of mutation towards deleterious mutants and s is the fitness cost of those deleterious mutants.
We can also compute the probability that no mutant is present at the start of vaccination:
When either
or r is small, p0 ≈ 1, and we can neglect the presence of preexisting mutants. Otherwise, we need to account for the possibility that one or more mutants are present at time t = 0, which we discuss in Section 5.5 below.
We now assume that there is zero mutant present at the start of vaccination. We are interested in the law of the first time tm at which a mutant appears. Because θU IUw(t) + θV IV w(t) is the flux of vaccine-escape mutants from the wild-type population, by the exponential formula for Poisson point processes, we have [27]:
We can numerically compute the probability density function fm of the first arrival time tm of a vaccine-escape mutant using
where Fm(t) is given by the auxiliary equation
with initial condition Fm(0) = 0, while we compute IUw(t) and IV w(t) using (12). The use of this auxiliary equation reduces computational effort by obtaining Fm(t) simultaneously with IUw(t) and IV w(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 tm ≥ 0 in a single infected host, that is, with density , (which is effectively zero as n becomes large). Then, (2) yields Im(tm) ≡ 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, , to the number of individuals infected with the mutant strain at time t after the mutant arrival time
. We approximate the rate of new infections,
by replacing the stochastic quantities
and Nn(t) by their deterministic approximations, giving the time-dependent birth rate
where SU (t), SV (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 dm = δ + αm + γ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 tm either goes extinct, or grows without bound (i.e. invades) with probability (see [24] and SI, Section 2)
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 tm using where
is obtained via the pair of auxiliary functions
and
[24] defined as follows:
and
, We then have
where
and we compute bm(t), SU (t), SV (t) and N (t), via (12). In practice, we cannot compute
; to obtain an approximation we evaluate
for sufficiently large t that
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 of a variant introduced at time t to characterize the distribution of the first time ti at which a mutant is introduced that successfully invades. By the thinning property of Poisson point processes, we have [27]:
where θU IUw(t) + θV IV w(t) is the flux of vaccine-escape mutants from the wild-type population. We compute numerically the probability density function gm of the first arrival time ti of a vaccine-escape mutant that successfully invades using
where Gm(t) is given by the auxiliary equation
with initial condition Gm(0) = 0 and computing IUw(t) and IV w(t) using (12).
Compare (22) with (17) and note that the probability that no vaccine-escape mutant will ever arise is
In contrast, the probability that no vaccine-escape mutant will ever invade is the larger probability
Note that
converges as t → ∞ to
which is nonzero, so that
that is, the probability of adaptation is 1 if and only if t 1→ (θU IUw(t) + θV IV w(t)) is not integrable. In other words, the probability of adaptation is 1 in the limit t → ∞ when the wild type is not driven to extinction by vaccination (i.e. ν < νc) which implies that there is an uninterrupted flux of mutation producing vaccine-escape variants. One of these mutants will eventually escape extinction and invade. Yet, the time needed for a successful variant to appear may be very long and we focus in the main text on
the probability of adaptation before time t (equation (10) and Figure 6).
5.4 Fixation of the variant (step 3)
Suppose now that the mutant strain successfully invades; we next consider the probability Pfix 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 Tm and Tw be the extinction times of mutant and wild-type strains, the probability of mutant fixation is thus ℙ{Tw < Tm} which we may decompose as
We again obtain estimates of
and
using the fact that conditional on SU (t), SV (t) and N (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
and the death rates are di = δ + αi + γ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 SV and SU. To quantify these epidemiological perturbations, we now approximate the density of susceptibles and total host density by the values of SU (t), SV (t) and N (t) obtained from the full deterministic system (1), which accounts for the presence of the mutant by replacing Im 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 tm. Then, at time tm, we use (1) where SU (tm), SV (tm), N (tm) and Iw(tm) are obtained from (12) and take (see Supplementary Information §4). Crucially, the initial density of the mutant depends on the probability of successful invasion of the mutant
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:
where
is the probability that an individual infected with strain i (i = m, w) present at tm 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,
.
In practice, we need two pairs of auxiliary equations to track the probability of extinction of both the wild type and the mutant
with
. 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).
Using (26) to approximate the probabilities
and
, we get
Differentiating this and taking
and
yields the following auxiliary equation for
:
with initial condition
. We estimate the fixation probability as
as above.
5.5 Invasion and Fixation with Standing Variation
We now come back to the probability of adaptation from standing variation at time t = 0, using the probability pk that there are k mutants present at time t = 0 and the estimates for the probabilities of invasion and fixation and
, of a single mutant arriving at time tm, taking tm = 0. Recall that pk is negative binomial with success probability
and
failures. Under the branching process approximation, the chain of infections started by each mutation will go extinct independently with probability
. The probability of invasion is then the probability that at least one line survives,
. Summing this over all possible values of k gives us the invasion probability from standing variation,
Recalling that the probability generating function for the number of mutants at time t = 0 is
which converges provided
, we see that the probability of invasion from standard variation is
Similarly, if there are k mutations at time t = 0, (26) gives us that , 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
. Substituting this for
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:
As previously, we obtain
by choosing t sufficiently large that
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 tm 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”.
Footnotes
E-mails:; AL - amaury.lambert{at}ens.fr; TD - tday{at}mast.queensu.ca; TLP - todd.parsons{at}upmc.fr;
We reworded several parts of the manuscript and Todd L Parsons was added as a coauthor