Advanced Search

See all authors and affiliations

Edited by Eugene Koonin, National Center for Biotechnology Information, National Library of Medicine, NIH, Bethesda, MD; received June 22, 2021; accepted November 18, 2021

The ongoing pandemic of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) raises an important question: who should we vaccinate first? Answering this question requires an analysis of both the short-term (epidemiological) and the long-term (evolutionary) consequences of targeted vaccination strategies. We analyze the speed of pathogen adaptation and the cumulative number of deaths in heterogeneous host populations to shed light on the effects of alternative vaccination strategies. This analysis shows that minimizing the speed of pathogen adaptation does not always minimize the number of deaths. This evaluation of both the epidemiological and evolutionary consequences of public health policies provides a practical tool to identify the best vaccination strategy.

The limited supply of vaccines against severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) raises the question of targeted vaccination. Many countries have opted to vaccinate older and more sensitive hosts first to minimize the disease burden. However, what are the evolutionary consequences of targeted vaccination? We clarify the consequences of different vaccination strategies through the analysis of the speed of viral adaptation measured as the rate of change of the frequency of a vaccine-adapted variant. We show that such a variant is expected to spread faster if vaccination targets individuals who are likely to be involved in a higher number of contacts. We also discuss the pros and cons of dose-sparing strategies. Because delaying the second dose increases the proportion of the population vaccinated with a single dose, this strategy can both speed up the spread of the vaccine-adapted variant and reduce the cumulative number of deaths. Hence, strategies that are most effective at slowing viral adaptation may not always be epidemiologically optimal. A careful assessment of both the epidemiological and evolutionary consequences of alternative vaccination strategies is required to determine which individuals should be vaccinated first.

The development of effective vaccines against severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) raises hope regarding the possibility of eventually halting the ongoing pandemic. However, vaccine supply shortages have sparked a debate about the optimal distribution of vaccination among different categories of individuals. Typically, infections with SARS-CoV-2 are far more deadly in older individuals than in younger ones (1). Prioritizing vaccination for older classes may thus provide a direct benefit in terms of mortality (2, 3). Yet, younger individuals are usually more active, and consequently, they may contribute more to the spread of the epidemic. Prioritizing vaccination for younger and more active individuals may thus provide an indirect benefit through a reduction of the epidemic size (4, 5). Earlier studies have compared alternative ways to deploy vaccination in heterogeneous host populations and showed that recommendation varies with the choice of the quantity one is trying to minimize (e.g., the cumulative number of deaths, the remaining life expectancy, or the number of infections) (3, 6, 7). The recommendation also varies with the properties of the pathogen and the efficacy of the vaccine (3, 4, 8). For SARS-CoV-2, the increase in mortality with age is such that the direct benefit associated with vaccinating more vulnerable individuals tends to overwhelm the indirect benefits obtained from vaccinating more active individuals (2, 3, 9, 10). However, some studies challenge this view and identified specific conditions where vaccinating younger and more active classes could be optimal (5, 7, 11, 12). A similar debate emerges over the possibility to delay the second vaccination dose to maximize the number of partially vaccinated individuals. A quantitative exploration of alternative vaccination strategies can help provide useful recommendations: a two-dose strategy is recommended when the level of protection obtained after the first dose is low and/or when vaccine supply is large (13⇓⇓–16).

Vaccine-driven evolution, however, could erode the benefit of vaccination and alter the above recommendations which are based solely on the analysis of epidemiological dynamics. Given that hosts differ both in their sensitivity to the disease and in their contribution to transmission, who should we vaccinate first if we want to minimize the spread of vaccine-adapted variants? The effect of alternative vaccination strategies on the speed of pathogen adaptation remains unclear. Previous studies of adaptation to vaccines focused on long-term evolutionary outcomes (17, 18). These analyses are not entirely relevant for the ongoing pandemic because what we want to understand first is the short-term consequence of different vaccination strategies (19). A few studies have discussed the possibility of SARS-CoV-2 adaptation following different targeted vaccination strategies but did not explicitly account for evolutionary dynamics (12, 20). A recent simulation study explored the effect of a combination of vaccination and social distancing strategies on the probability of vaccine-driven adaptation (21). This model, however, did not study the impact of targeted vaccination strategies on the speed of adaptation.

Here we develop a theoretical framework based on the analysis of the deterministic dynamics of multiple variants after they successfully managed to reach a density at which they are no longer affected by the action of demographic stochasticity. We study the impact of different vaccination strategies on the rate of change of the frequency of a novel variant, which allows us to quantify the speed of virus adaptation to vaccines. Numerical simulations tailored to the epidemiology of SARS-CoV-2 confirm the validity of our approximation of the strength of selection for vaccine-adapted variants.

We are interested in tracking the frequency *p _{m}* of hosts infected by the vaccine-adapted variant among all the infected hosts. It is possible to show that under a broad range of conditions one can approximate the dynamics of the vaccine-adapted variant frequency asp˙m≈pm(1−pm)S(t),[1]where S(t) is the selection coefficient on the vaccine escape mutation. This selection coefficient measures the rate of change of the logit of the frequency of the vaccine-adapted variant [i.e., ln (pm/(1−pm))] and provides a relevant measure for the speed at which the viral population is adapting (

*Materials and Methods*).

Targeted vaccination strategies aim to preferentially vaccinate hosts according to specific epidemiological characteristics. For instance, we could target hosts that have more contacts or are more at risk for a severe disease. In our model, we therefore introduce some heterogeneity among hosts. As a result, from the point of view of the parasite, the quality of the host may differ among infected hosts, and this variation is likely to affect the dynamics of vaccine-adapted variants. To quantify host quality, we use the concept of reproductive value, a key concept in demography and evolutionary biology (22⇓–24). Reproductive value measures how much a virus infecting a given class of hosts will contribute to the future of the viral population. Our general mathematical analysis allows us to take the difference in host quality into account when calculating the selection coefficient S(t) (

*Materials and Methods*).

We use this approach to analyze the speed of adaptation during the ongoing pandemic of SARS-CoV-2 under different scenarios. We use an epidemiological model tailored to the biology of SARS-CoV-2 (

*Materials and Methods*). However, it is important to keep in mind that due to simplifying assumptions and uncertainty about parameter values, our results cannot be translated directly into public health recommendations without further investigations (

*Discussion*). Nonetheless, our theoretical framework gives clear foundations for future applied work and captures some of the most salient features of the COVID-19 pandemic. In particular, we introduce a time-varying parameter

*c*(

*t*), which measures the intensity of nonpharmaceutical interventions (NPI). We assume that the epidemic is initially controlled by NPI, which yields successive epidemic waves before the deployment of vaccination at

*t*= 150 d. We use this model to explore the effect of two different forms of heterogeneity on the speed of SARS-CoV-2 adaptation.

In the first scenario we assume that hosts differ in their ability to mix and thus to transmit the disease. More specifically, following the model used by ref. 12, we assume that some hosts (

*L*) have a low number of social interactions, while other hosts (

*H*) have a higher number of contacts. These two types of hosts can be thought as corresponding to the older and younger halves of the population. The increased rate of social interactions among

*H*hosts is captured by a parameter M≥1. Susceptible hosts are initially naive (

*S*and

^{L}*S*), but they can become vaccinated (S^L and S^H) at rates

^{H}*ν*and

^{L}*ν*, respectively. When vaccinated, hosts have a lower probability to become infected (rσ measures the efficacy at blocking infection), and if they become infected, they have a lower probability to transmit the virus (rτ measures the efficacy at blocking transmission) and to die from the infection (rμ measures the efficacy at reducing mortality). Viral adaptation, however, can erode these benefits. We consider different viral strains characterized by an escape trait

^{H}*e*which takes values between 0 (no escape) and 1 (full escape). The capacity of a variant to reduce the effect of the vaccine on transmissibility (infectivity) is captured by a function Eτ(e) for transmissibility (and Eσ(e) for infectivity), which allows us to quantify the overall ability of the virus to escape the protective effects of the vaccine as E(e)=Eτ(e)Eσ(e). Note that the capacity of a variant to reduce mortality does not affect the strength of selection in our model (i.e., the duration of infection is affected neither by the variant nor by the vaccine).

In

*Materials and Methods*, we derive a simple approximation for the strength of selection acting on the vaccine-adapted variant:S(t)=(1−c(t)) β ΔE (S^L+M2S^H),[2]where ΔE refers to the change in vaccine escape ability caused by the mutation. This tells us that the intensity of selection depends on 1) the ability of the virus mutant to escape the protective effects of vaccine, 2) the densities of uninfected hosts (both

*L*and

*H*) who have been vaccinated, and 3) the relative number of contacts of each class of hosts. Note that the epidemiological impact of a higher contact rate (M>1) translates into a magnified selective impact (M2). Thus, if we have to choose between vaccinating

*L*and

*H*hosts, targeting

*H*hosts is expected to select more strongly for the vaccine-adapted variant. Fig. 1

*B*confirms that this approximation captures very well the temporal dynamics of the vaccine-adapted variant. In particular, the simulations confirm that targeted vaccination of the

*L*hosts slows down the rate of adaptation of the virus.

Scenario 1: heterogeneity in contact numbers and vulnerability. (

*A*) A graphical presentation of the epidemiological life cycle with

*L*hosts who are more vulnerable to the disease and

*H*hosts who have a higher number of contacts. Infected hosts are indicated with a light red shading, and vaccination is indicated with a bold circle border. The force of infection on naive hosts is noted Λi=hi+h^i (see

*Materials and Methods*and Table 1 for additional details on this model). (

*B*) Dynamical change of the frequency

*p*of the vaccine-adapted mutant for two distinct targeted vaccination strategies: 1) mostly

_{m}*L*hosts are vaccinated (blue lines) and 2) mostly

*H*hosts are vaccinated (red lines). The full lines indicate the exact numerical computation, and the dashed lines indicate the approximation obtained from Eq.

**2**. The gray areas indicate the period where NPI were used to control the epidemic (c(t)=0.7 with NPI). (

*C*) Incidence of the epidemic (fraction of the total host population that is infected) in the absence of vaccination (dotted black line) or under the two alternative vaccination strategies used in

*B*(blue and red lines). (

*D*) Cumulative number of deaths (fraction of the total host population) in the absence of vaccination (dotted black line) or under the two alternative vaccination strategies used in

*B*(blue and red lines).

Of course, the choice of the vaccination strategy should not be based solely on the reduction of the speed of adaptation to vaccines. Indeed, the best way to limit the spread of vaccine escape mutations would be to adopt the worst epidemiological strategy: avoiding the use of vaccines. Yet, we urgently need vaccines to save lives and halt the current pandemic. We can use our numerical simulations to study the consequences of distinct targeted vaccination strategies on the total number of cases and on mortality (

*Materials and Methods*). Fig. 1

*C*shows that targeting

*L*hosts is expected to increase the number of cases because

*H*hosts contribute more to the spread of the disease. Yet, Fig. 1

*D*shows that targeting

*L*hosts is expected to decrease the cumulative number of deaths after some time because

*L*hosts (i.e., older individuals) are also associated with higher risks of dying from the infection. Hence, targeting

*L*hosts makes sense both for epidemiological and evolutionary reasons.

We explored the robustness of the above results for a range of alternative scenarios. First, we note that as expected from our analytic approximation, the use of a transmission-blocking vaccine (instead of an infection-blocking vaccine) yields very similar outcomes (compare Fig. 1 and

*SI Appendix*, Fig. S1). Second, we show in

*SI Appendix*, Fig. S2 that evolution amplifies the increase in the cumulative number of deaths when

*H*hosts are vaccinated compared to a scenario without viral adaptation. Indeed, the spread of a vaccine-adapted variant drives a large epidemic wave in vaccinated populations. This evolutionary effect is maximized for intermediate values of the speed of the vaccination rollout because when vaccine rollout is very fast, the vaccine-adapted variant is rapidly favored, whatever the targeted vaccination strategy (

*SI Appendix*, Fig. S3). Finally, we note that maintaining social distancing for longer can substantially decrease the speed of adaptation (

*SI Appendix*, Fig. S4).

In our second scenario, we assume that the heterogeneity among hosts is determined by differential strength of immunity induced by distinct vaccination status. We distinguish between unvaccinated hosts (

*S*), hosts partially vaccinated with one dose (S^I), and hosts “fully vaccinated” with two doses (S^II). Using the same approach as before, we obtain the following expression for the strength of selection acting on the strength of selection on the vaccine-adapted variant:S(t)=(1−c(t)) β (ΔEI S^I+ ΔEII S^II).[3]

Eq.

**3**is very similar to Eq.

**2**, but now we have to account for the fact that the escape mutation has different effects in each class. Hence, the influences of an increase in the densities of hosts vaccinated by a single or two doses of vaccines are weighted by ΔEI and ΔEII, respectively. A single vaccine dose is likely to induce a lower protection against the virus (i.e., EI>EII), but this does not necessarily imply that ΔEI>ΔEII. In fact, we can show that if the vaccine is acting on a single step of the virus’ life cycle (e.g., only blocking infection), we expect ΔEII>ΔEI. Delaying the acquisition of the second dose will have two effects: 1) a lower density S^II of fully vaccinated hosts decreases the more intense selection imposed by these hosts but 2) delaying the second dose allows for more hosts to be vaccinated, and the increase in S^I may result in stronger selection for the vaccine-adapted variants. We show in Fig. 2

*B*that this second effect can be more important than the first one, and delaying the second dose can result in faster adaptation. However, Fig. 2

*D*shows that delaying the second dose may reduce the cumulative number of deaths because a larger fraction of the population would benefit from the protection of the vaccine (but higher rates of vaccination rollout can reverse this effect on mortality;

*SI Appendix*, Figs. S5 and S6). Hence, in contrast to the previous scenario, the strategy that maximizes the speed of adaptation may result in a lower mortality. The contrast between our two scenarios illustrates the necessity to quantify both the epidemiological and the evolutionary consequences of different targeted vaccination strategies to identify the optimal way to distribute vaccines.

Scenario 2: heterogeneity in the number of vaccination doses. (

*A*) A graphical presentation of the epidemiological life cycle where the superscripts I and II refer to the first and second doses of vaccine. Infected hosts are indicated with a light red shading, and vaccination is indicated with a bold circle border. The force of infection on naive hosts is noted Λi=hi+h^i (see

*Materials and Methods*and Table 2 for additional details on this model). (

*B*) Dynamical change of the frequency

*p*of the vaccine-adapted mutant for two distinct targeted vaccination strategies: 1) vaccinated hosts receive two doses sequentially (purple lines) and 2) a single dose is used for each host (orange lines). The full lines indicate the exact numerical computation, and the dashed lines indicate the approximation obtained from Eq.

_{m}**3**. The gray areas indicate the period where NPI were used to control the epidemic (c(t)=0.7 with NPI). (

*C*) Incidence of the epidemic (fraction of the total host population that is infected) in the absence of vaccination (dotted black line) or under the two alternative vaccination strategies used in

*B*(purple and orange lines). (

*D*) Cumulative number of deaths (fraction of the total host population) in the absence of vaccination (dotted black line) or under the two alternative vaccination strategies used in

*B*(purple and orange lines).

The speed of the spread of SARS-CoV-2 variants has baffled the scientific community (25, 26). Despite a relatively small mutation rate (27, 28), SARS-CoV-2 has the ability to produce mutations with variable phenotypic effects that fuel the adaptation to human populations. The growing concern regarding the ability of the virus to escape host immunity calls for tools allowing us to anticipate the speed of the spread of vaccine escape mutants. We show here that heterogeneity in the behavior (scenario 1) and/or immune status (scenario 2) of the human host can induce variation in the strength of selection for vaccine escape mutations in the virus population. We contend that it is important to quantify this variation because it could be used to carry out targeted vaccination strategies that, for a given vaccination coverage, could limit the speed of adaptation of the virus.

We show that targeted vaccination on older hosts which are associated with lower number of contacts but higher risks of mortality may be a good strategy to reduce both the spread of the vaccine escape variant and the cumulative number of deaths. Ref. 12 used a different approach to identify vaccination strategies that could reduce what they call “vaccine escape pressure,” a quantity proportional to the density of infected hosts who are vaccinated. In contrast, we show that the strength of selection on the vaccine-adapted mutant is proportional to the density of susceptible hosts who are vaccinated. Their analysis relies on the assumption that the incidence of the infection remains small (i.e., no depletion of susceptible hosts), and they do not track explicitly the rate of spread of a vaccine-adapted variant. They conclude that vaccinating most of the vulnerable hosts and few of the mixers could be the most risky for vaccine escape. Yet, it is difficult to evaluate how the vaccine escape pressure criteria used in ref. 12 may affect the speed of viral adaptation. A high incidence among vaccinated hosts may speed up viral adaptation because a vaccine-adapted variant generated by mutation is more likely to escape extinction in a vaccinated host. However, a quantification of this effect would require an explicit description of the interplay between within-host selection and demographic stochasticity. In other words, their model focuses on the process that limits the emergence of vaccine-adapted variants, while our analysis focuses on the strength of selection after emergence. It would thus be particularly interesting to explore the robustness of our results with a model that would account for the effects of demographic stochasticity and within-host selection on the emergence of new variants.

We also discuss the effect of delaying the second dose of the vaccine on viral adaptation and on mortality. In a recent model, ref. 29 found that imperfect immunity induced by a single dose may lead to stronger within-host selection for vaccine escape variants. This is the same argument used by ref. 12, where the infection of imperfectly immunized hosts may speed up viral adaptation. However, as discussed above, this effect relies on the interplay between demographic stochasticity and within-host selection. In contrast, we focus on between-host selection and ask whether vaccine-adapted variants can increase in frequency at the population level. We contend that once vaccine-adapted variants reach a significant fraction of the population, the fate of those mutations will be driven by between-host selection. Our analysis clarifies the balance between the effects of mutations in different types of hosts (i.e., the relative magnitudes of ΔEI and ΔEII) and the quantity of the different types of hosts (i.e., the relative densities of hosts with one or two doses of the vaccine). We show that a higher speed of adaptation may be the price to pay for a reduced number of deaths (Fig. 2

*B*). Indeed, delaying the second dose allows for protecting (albeit partially) a larger fraction of the population (see ref. 15 for an exploration of this effect). This positive effect can outweigh the negative consequences of an erosion of vaccine efficacy due to viral adaptation.

Interestingly, we found that combining vaccination and NPI can delay the rise of vaccine-adapted variants (

*SI Appendix*, Figs. S4 and S7). Our analysis clarifies the origin of the positive effect of NPI discussed in previous studies (21, 30). In our deterministic model, this effect emerges from the reduction in the strength of between-host selection due to lower opportunities of transmission (i.e., see the effect of larger values of

*c*(

*t*) in Eqs.

**2**and

**3**).

Vaccination is urgently needed to control the SARS-CoV-2 pandemic, but the limited supply of vaccines is raising major ethical and practical issues. Public health policies need to strike a balance between social, ethical, and short-term epidemiological considerations. Our work illustrates that the long-term evolutionary consequences of specific vaccination strategies also need to be considered and evaluated using quantitative models. Indeed, viral adaptation could erode the efficacy of vaccines, and targeted vaccination may provide a way to delay this adaptation. Yet, as illustrated with the second scenario, a strategy that minimizes the cumulative number of deaths may not necessarily minimize the speed of adaptation. Hence, as for any therapeutic interventions that may result in the evolution of pathogen resistance, the identification of an optimal vaccination strategy that reduces the death toll of the pandemic requires specific models accounting for both the epidemiology and the evolution of the virus (31⇓–33). Because our model relies on several simplifying assumptions and because our knowledge of the biology of SARS-Cov-2 and of several key parameter values (e.g., vaccine efficacy, virulence, and contact rates) remains imperfect, our model cannot be used directly to make quantitative public health recommendations. Nonetheless, our framework lays a clear conceptual foundation to analyze the consequences of targeted vaccination strategies. In order to make more precise applied predictions, it would be interesting to investigate how other realistic factors (such as age structure, difference in transmissibility among host classes, or alternative vaccination schedules) may affect our results. Furthermore, the present work could be readily extended to combine the two forms of heterogeneities in the same model to allow for alternative ways to distribute the vaccines (e.g., two doses for

*L*individuals and one dose for

*H*individuals). It would also be possible to use the same framework to account for other factors that have been shown to affect the outcome of vaccination strategy like assortative mixing (4) and compensatory behavior after vaccination (34). In all these scenarios, our framework could be used to identify which strategy manages to strike the right balance between the epidemiological and the evolutionary consequences of targeted vaccination strategies.

We first give a general overview of the method used to calculate the selection coefficient in structured host populations. The dynamics of hosts infected by pathogen strain

*i*can be captured by a matrix Ri collecting the transition rates between host classes. Assuming that mutations have small phenotypic effects (i.e., em=ew+ε), we can write the change in frequency of the mutant strain asdpmdt=εpm(1−pm) v⊤dRmdemf+O(ε2),[4]where v⊤ is the vector of reproductive values and

**f**is the vector of class frequencies. These vectors are conormalized such that v⊤f=1 and satisfy the following dynamical equations:dfdt=Rwf−λ(t)f,[5]dv⊤dt=−v⊤Rw+λ(t)v⊤,[6]

where Rw is the transition matrix for the wild-type strain and λ(t) is the per-capita growth rate of the resident population at time

*t*(see refs. 24, 35 for a more detailed description). The dRm/dem term in Eq.

**6**refers to the differentiation of each element of the transition matrix Rm with respect to the effect of the mutant. For small

*ε*, the mutant frequency

*p*changes slowly compared to the ecological variables

_{m}**f**and

**v**, and we can use a quasi-equilibrium approximations obtained by setting the right-hand sides of Eqs.

**5**and

**6**to zero. This allows us to obtain analytical expressions for the class frequencies and reproductive values and thus to calculate the selection coefficient for a specific life cycle (scenario 1 vs. scenario 2). Note that although the weak selection assumption (small

*ε*) is driving the separation of time scale, the approximation remains good when selection is strong as discussed in the two scenarios below and shown in Figs. 1 and 2. The Mathematica notebooks used to generate Figs. 1 and 2 and

*SI Appendix*, Figs. S1–S7 are accessible from the following data repository, https://doi.org/10.5281/zenodo.5803062.

We assume that susceptible hosts in class

*k*(where

*k*=

*L*or

*H*) are vaccinated at rate

*ν*representing the speed of vaccination rollout in that host class. We note M>1 as the relative number of contacts of

^{k}*H*hosts compared to

*L*hosts and ρτ (ρσ) as the relative transmissibility (susceptibility) of vaccinated hosts compared to naive hosts of the same class. Both ρτ and ρσ are functions of the vaccine escape trait. With these assumptions, the force of infection of a pathogen strain

*i*due to naive infected hosts is hi=βiLIiL+MβiHIiH, and h^i=ρτ(ei)(βiLI^iL+MβiHI^iH) for vaccinated infected hosts. Note that vaccinated hosts are indicated by a “hat” (denoting protection). Hosts in class

*k*infected by pathogen strain

*i*eventually leave the class at rate δik (δ^ik for vaccinated hosts) and can either recover or die. We assume that the probability μik (μ^ik) of dying after leaving the class Iik (I^ik) may depend on the host class

*k*and pathogen strain

*i*. We track the cumulative number of deaths

*D*. This quantity can be used to compare the efficacy of different vaccination strategies. Note that the probabilities μik and μ^ik have no impact on evolutionary dynamics because these events occur when the host is assumed to be no longer infectious, and consequently, they do not affect pathogen fitness.

This yields the following dynamical system (see also Fig. 1

*A*):S˙L=−νLSL−(hi+h^i)SLS^˙L=νLSL−(hi+h^i)ρσ(ei)S^LS˙H=−νHSH−M(hi+h^i)SHS^˙H=νHSH−M(hi+h^i)ρσ(ei)S^HI˙iL=(hi+h^i)SL−δiLIiLI^˙iL=(hi+h^i)ρσ(ei)S^L−δ^iLI^iLI˙iH=M(hi+h^i)SH−δiHIiHI^˙iH=M(hi+h^i)ρσ(ei)S^iH−δ^iHI^iHD˙=∑i(μiLδiLIiL+μ^iLδ^iLI^iL+μiHδiHIiH+μ^iHδ^iHI^iH).[7]

We analyze this general model under two simplifying but reasonable assumptions.

1. We assume that the pathogen strains only differ through their effect on the parameters ρτ and ρσ (that is, we only look at vaccine escape mutations, not mutations that can also affect transmissibility or virulence).

2. We assume that host classes

*L*and

*H*only differ through their number of contacts, so that δiL=δ^iL=δiH=δ^iH=δ and βiL=βiH=β(1−c(t)), where

*β*is the baseline transmissibility and 1−c(t) captures the effect of NPI aimed at controlling the epidemic by reducing transmission. The parameter

*c*(

*t*) varies between 0 and 1 and quantifies the intensity of the control, which may vary over time as observed during the COVID-19 pandemic.

The latter assumption implies that the duration of infection is the same in all classes, but the effect of vaccination on the mortality of the different classes of hosts can be captured through the probabilities

*μ*, μ^L,

^{L}*μ*, and μ^H (again, note that there is no influence of the pathogen genotype on disease outcome). For instance, in our simulations, we assume that

^{H}*L*hosts tend to have fewer contacts but a higher mortality risk, while

*H*hosts have more contacts but a lower mortality risk. This may reflect the observed differences between age classes.

The transition matrix Ri is the 4 × 4 matrix of per-capita transition rates of the pathogen between the four different types of hosts, given by Eq.

**8**.Ri=(βLSL−δLβLρτ(ei)SLβLρσ(ei)S^LβLρτ(ei)ρσ(ei)S^L−δ^LMβLSHMβLρτ(ei)SHMβLρσ(ei)S^HMβLρτ(ei)ρσ(ei)S^H MβHSLMβHρτ(ei)SLMβHρσ(ei)S^IMβHρτ(ei)ρσ(ei)S^LM2βHSH−δHM2βHρτ(ei)SHM2βHρσ(ei)S^HM2βHρτ(ei)ρσ(ei)S^H−δ^H).[8]

We are interested in the dynamics of the frequency of the vaccine escape mutant, which ispm=ImL+ImH+I^mL+I^mH∑i(IiL+IiH+I^iL+I^iH).[9]

The dynamics of

*p*can be calculated by plugging the expressions of Rw and Rm into Eqs.

_{m}**4**–

**6**. After some rearrangements, we obtainS(t)=dρσdemhw(v^LS^L+Mv^HS^H) +dρτdemh^w(vLSL+ρσ(ew)v^LS^L+MvHSH +Mρσ(ew)v^HS^H),[10]where the vector vT=(vLv^LvHv^H) collects the reproductive values of an individual resident pathogen in classes

*I*, I^L,

^{L}*I*, and I^H, respectively. Note that this result only depends on assumption 1 above.

^{H}It is possible to simplify the expression of the selection coefficient by treating the reproductive values as fast variables. In particular, using our assumption 2, this leads to the following quasi-equilibrium approximations:v^LvL=v^HvH=ρτ(ew), vHvL=M, v^HvL=Mρτ(ew).[11]

Similarly, we have the following quasi-equilibrium approximations for the class frequencies, which give the fraction of infected individuals in a given class:f^LfL=ρσ(ew)S^LSL, fHfL=MSHSL, f^HfL=Mρσ(ew)S^HSL.[12]

Together with the normalization condition vLfL+v^Lf^L+vHfH+v^Hf^H=1, we can use these relationships to obtainS(t)=εβ(1−c(t))(S^L+M2S^H)d(ρσρτ)dem∣em=ew.[13]

To recover Eq.

**2**, we use the notations Eτ=ρτ, Eσ=ρσ, E=EτEσ, andΔE=εdEdem∣em=ew,which is the first-order approximation of the difference E(em)−E(ew).

In our applications, we use a linear model of vaccine escape:ρτ(ei)=1−rτ(1−ei),[14]ρσ(ei)=1−rσ(1−ei),[15]where rτ and rσ give the vaccine efficacy in the absence of vaccine escape mutation (i.e.,

*e*= 0). When

_{i}*e*= 1 (full vaccine escape), the vaccine offers no reduction in transmissibility and susceptibility (ρσ=ρτ=1). As explained above, we assume that

_{i}*L*hosts have a higher risk of mortality due to the disease and note D as the relative increase in mortality of

*L*hosts vs.

*H*hosts,

*μ*as the baseline mortality probability, and ρμ=1−rμ as the reduction in mortality due to the vaccine (which we assume independent of host classes and of pathogen genotype). We thus haveμL=D μμ^L=ρμ D μμH=μμ^H=ρμ μ,

with D>1 and 0<ρμ<1.

Initial conditions used in Fig. 1

*B–D*are as follows: SL(0)=SH(0)=1/2, S^L=S^H=0, IL(0)=IH(0)=10−6, I^L(0)=I^H(0)=10−6, D(0)=0, and pm(0)=10−3. The intensity of NPI varies with time (c(t)=0.7 when t∈[40,140] and t∈[150,250],

*c*(

*t*) = 0 otherwise). Vaccination starts at

*t*= 150, and the other parameters used in Fig. 1

*B–D*are listed in Table 1.

Main parameters and default values used in scenario 1

For our second scenario, we consider three classes of susceptible hosts: unvaccinated (

*S*), vaccinated with one dose (SI), and vaccinated with two doses (SII). Unvaccinated susceptible hosts can be given a first dose of vaccine at rate νI. Susceptible hosts that have received one dose can be given a second dose at rate νII. With one dose, the relative transmissibility (susceptibility) of vaccinated hosts with respect to pathogen strain

*e*is ρτI(ei) (ρσI(ei)). With two doses, we use the notation ρτII(ei) and ρσII(ei). Apart from these assumptions, the life cycle is similar to the one used for scenario 1, and we have the following dynamics (see also Fig. 2

_{i}*A*):S˙=−νIS−(hi+h^i)SS^˙I=νIS−νIIS^I−(hi+h^i)ρσI(ei)S^IS^˙II=νIIS^I−(hi+h^i)ρσII(ei)S^III˙i=(hi+h^i)S−δiIiI^˙iI=(hi+h^i)ρσI(ei)S^I−δ^iII^iII^˙iII=(hi+h^i)ρσII(ei)S^II−δ^iIII^iIID˙=∑i(ciδiIi+c^iIδ^iII^iI+c^iIIδ^iIII^iII)),[16]where the forces of infection by virus strain

*i*are hi=βiIi and h^i=ρτI(ei)βiII^iI+ρτII(ei)βiIII^iII. For simplicity, we will also assume, as in scenario 1, that βi=βiI=βiII=β(1−c(t)) and δi=δ^iI=δ^iII=δ, so that 1) hosts only differ through the parameters ρτ and ρσ and 2) the viral strains only differ through the parameters ρτ and ρσ. We also assume that μi=μ, μ^iI=μ^I, and μ^iII=μ^II to account for potential differences between mortality rates between different classes of hosts (but no influence of the pathogen genotype).

With these assumptions, the matrix Ri isRi=β(1−c(t))(SρτI(ei)SρτII(ei)SρσI(ei)S^IρτI(ei)ρσI(ei)S^IρτII(ei)ρσI(ei)S^IρσII(ei)S^IIρτI(ei)ρσII(ei)S^IIρτII(ei)ρσII(ei)S^II)−δJ,[17]where

**J**is the 3 × 3 identity matrix, and the quasi-equilibrium relationships for class frequencies and reproductive values arevIv0=ρτI(ew), vIIv0=ρτII(ew)

andfIf0=ρσI(ew)S^IS, fIIf0=ρσII(ew)S^IIS,

where v=(v0vIvII) and f=(f0fIfII). Together with the normalization condition v0f0+vIfI+vIIfII=1, these relationships allow us to rearrange Eq.

**4**to obtainS(t)=εβ(1−c(t))(S^Id(ρσIρτI)dem∣em=ew+S^IId(ρσIIρτII)dem∣em=ew),[18]

which is Eq.

**3**using the same notations as in scenario 1.

We use the same linear model of vaccine escape as in scenario 1, but we allow for different vaccine efficacies depending on the number of doses:ρτI(ei)=1−rτI(1−ei),[19]ρσI(ei)=1−rσI(1−ei),[20]ρτII(ei)=1−rτII(1−ei),[21]ρσII(ei)=1−rσII(1−ei),[22]where rτ and rσ give the vaccine efficacy in the absence of vaccine escape mutation (i.e.,

*e*= 0). When

_{i}*e*= 1 (full vaccine escape), the vaccine offers no reduction in transmissibility and susceptibility (ρσ=ρτ=1). We assume that vaccination can also protect against disease-induced mortality, and we define rμI and rμII so that μ^iI=μ(1−rμI) and μ^iII=μ(1−rμII).

_{i}Initial conditions used in Fig. 2

*B–D*are as follows: S(0)=1, S^I(0)=S^II(0)=0, I(0)=I^I(0)=I^II(0)=10−6, D(0)=0, and pm(0)=10−3. The intensity of NPI varies with time (c(t)=0.7 when t∈[40,140] and t∈[150,250],

*c*(

*t*) = 0 otherwise). Vaccination starts at

*t*= 150, and the other parameters used in Fig. 2

*B–D*are listed in Table 2.

Main parameters and default values used in scenario 2

Mathematica notebooks have been deposited in Zenodo https://doi.org/10.5281/zenodo.5803062 (36).

Several countries recommend a third dose of vaccination. Our scenario 2 could be readily modified to account for this new class of vaccinated individuals and this would yield a third term, ΔEIII S^III, in Eq.

**3**.

We thank the Editor and two anonymous reviewers for helpful and constructive comments. This work was funded by grants ANR-16-CE35-0012 “STEEP” to S.L. and ANR-17-CE35-0012 “EVOMALWILD” to S.G. from the Agence Nationale de la Recherche.

↵

^{1}S.G. and S.L. contributed equally to this work.

Author contributions: S.G. and S.L. designed research, performed research, and wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2110666119/-/DCSupplemental.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Thank you for your interest in spreading the word on PNAS.

NOTE: We only request your email address so that the person you are recommending the page to knows that you wanted them to see it, and that it is not junk mail. We do not capture any email address.

Submit

Feedback Privacy/Legal

Copyright © 2022 National Academy of Sciences. Online ISSN 1091-6490. PNAS is a partner of CHORUS, CLOCKSS, COPE, CrossRef, ORCID, and Research4Life.