The speed of vaccination rollout and the risk of pathogen adaptation

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, I_i := I_Ui + I_{V i}, which yields: 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: α_i, β_i, γ_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: 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 (e_wR_w > 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 e_w and ω) yield lower threshold values for the speed of vaccination. Imperfect vaccines (i.e., higher values of e_w 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 f_m(t) of the arrival time t = t_m of the vaccine-escape mutant is approximated by: 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 ν 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.

• Download figure
• Open in new tab

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: ν = 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, e_w = 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 d_m = δ + α_m + γ_m, introduced at time t_m, successfully invades (see [24] and SI, Section 2) is: 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)).

• Download figure
• Open in new tab

Figure 4: Probability of invasion of the vaccine-escape mutant increases with t_m.

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: e_m = 1, β_m = 10, α_m = 1. The fast mutant: e_m = 1, β_m = 30, α_m = 5.02. The probability of invasion in the limit t_m → ∞ (see equation (9)) is indicated with the dashed black line. Other parameter values: λ = 0.01, δ = 0.01, γ = 1, β_w = 20, α_w = 1, e_w = 0.03, ω = 0.05, p = 0.1.

Taking t_m → ∞ 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 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 , 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 r_i, not the per-generation reproduction ratio R_i [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 t_m 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).

• Download figure
• Open in new tab

Figure 5: Probability of fixation of the vaccine-escape mutant may be low when t_m 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: ν = 0.1 (top), 0.2 (middle) and 0.3 (bottom). The slow mutant: e_m = 1, β_m = 10, α_m = 1. The fast mutant: e_m = 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, e_w = 0.03, ω = 0.05, p = 0.1, n = 10⁸.

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):

• Download figure
• Open in new tab

Figure 6: Probability of adaptation and the speed of vaccination rollout.

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, e_w = 0.03, β_m = 10, α_m = 1, e_m = 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).

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 I_w = I_Uw + I_{V w}, we get from (12): A new endemic equilibrium will thus be approached after vaccination if and only if the growth rate r_w = β_w(S_U + e_wS_V)/N − (δ + α_w + γ_w) is positive, or equivalently if the effective per-generation reproduction ratio is larger than 1, when S_U and S_V are taking their stationary values and in the absence of infection. Computing these values (see SI, Section 1) shows that if and only if e_wR_w > 1 or ν > ν_c, where Thus we see that if e_wR_w < 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 e_wR_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 θ_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 d_m = δ + α_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, p₀ ≈ 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 t_m at which a mutant appears. Because θ_U I_Uw(t) + θ_V I_{V 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 f_m of the first arrival time t_m of a vaccine-escape mutant using where F_m(t) is given by the auxiliary equation with initial condition F_m(0) = 0, while we compute I_Uw(t) and I_{V w}(t) using (12). The use of this auxiliary equation reduces computational effort by obtaining F_m(t) simultaneously with I_Uw(t) and I_{V 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 t_m ≥ 0 in a single infected host, that is, with density , (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, , 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 Nⁿ(t) by their deterministic approximations, giving the time-dependent birth rate 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 = δ + α_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 t_m 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 t_m using where is obtained via the pair of auxiliary functions and [24] defined as follows: and , We then have where and we compute b_m(t), S_U (t), S_V (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 t_i at which a mutant is introduced that successfully invades. By the thinning property of Poisson point processes, we have [27]: where θ_U I_Uw(t) + θ_V I_{V w}(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-escape mutant that successfully invades using where G_m(t) is given by the auxiliary equation with initial condition G_m(0) = 0 and computing I_Uw(t) and I_{V 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 I_Uw(t) + θ_V I_{V 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 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 ℙ{T_w < T_m} which we may decompose as We again obtain estimates of and using the fact that conditional on S_U (t), S_V (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 d_i = δ + α_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 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 (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 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, .

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 p_k 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 t_m, taking t_m = 0. Recall that p_k 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 ⁰(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 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.