Abstract
Several vaccines with different efficacies and effectivenesses are currently being distributed across the world to control the COVID19 pandemic. Having enough doses from the most efficient vaccines in a short time is not possible for all countries. Hence, policymakers may propose using various combinations of available vaccines to control the pandemic with vaccineinduced herd immunity by vaccinating a fraction of the population. The classic vaccineinduced herdimmunity threshold suggests that we can stop spreading the disease by vaccinating a fraction of the population. However, that classic threshold is defined only for a single vaccine and may be invalid and biased when we have multivaccine strategies for a disease or multiple variants, potentially leading policymakers to suboptimal vaccineallocation policies. Here, we determine which combination of multiple vaccines may lead to herd immunity. We show that simplifying the problem and considering the vaccination of the population as a singlevaccine strategy whose effectiveness is the sample mean of all effectivenesses would not be ideal, because many multivaccine strategies with a smaller herdimmunity threshold can be proposed. We show that the herdimmunity threshold may vary due to changes in vaccineuptake proportions. Moreover, we propose methods to determine the optimal combination of multiple vaccines in order to achieve herd immunity and apply our results to the issue of multiple variants. In addition, we determine a condition for reaching herd immunity in the presence of new emerging variants of concern. We show by example that new variants could influence our estimation of the vaccination reproduction number. It follows that the herdimmunity threshold must be updated not only when multivaccine strategies are used but also when multiple variants coexist in the population.
Introduction
Following the rapid worldwide spread of the coronavirus SARSCoV2, which first emerged in late December 2019 in Wuhan, China, a wide range of policies and interventions have been followed by countries in order to respond to the pandemic^{1,2}. Multiple waves of the pandemic appeared all around the world, and multiple variants have emerged. Vaccination is one of the main policies to control this worldwide pandemic^{3}. However, the vaccine stockpile is not sufficient to vaccinate the entire global population with any single vaccine, so using multiple vaccines is be the best chance to control the spread of the epidemic.
Herd immunity is a crucial determinant for publichealth policymakers to contain and potentially eradicate an infectious disease. It plays an important role in controlling the pandemic once a sufficiently high proportion of the population gains immunity through vaccination or infection. The herdimmunity threshold is the proportion of a population that must be vaccinated to stop an infectious disease from spreading. It can be approximated using the basic reproduction number, \(\text {R}_0\), which is an epidemic measure used to describe the contagiousness of infectious diseases. It is defined as average number of new infections caused by a single infectious individual in a completely susceptible population^{4,5}. Although \(\text {R}_0\) is affected by various biological, environmental, sociological and behavioural factors, it is generally reported as a single numeric value with a straightforward interpretation that incidence of the infectious disease is expected to decline if \(\text {R}_0\) is less than one and to increase if \(\text {R}_0\) is larger than one^{5,6,7}.
Suppose that a specific infectious disease is distributed homogeneously in the population and a single vaccine provides protection against it. Let p be the fraction of the susceptible individuals who get vaccinated, with vaccinees selected randomly with the same probability among susceptible individuals. We assume the vaccine provides lifelong immunity with no reinfection and protects only a proportion E among those vaccinated. Then the vaccination reproduction number is
where \(E\in (0,1]\) represents the vaccine effectiveness^{8,9}. If the vaccine is perfect and provides \(100\%\) immunity, then its effectiveness is \(E=1\), but if it is imperfect and provides only partial immunity, then \(E<1\). The population reaches herd immunity, with incidence of the infection declining, if \(\text {R}_V\le 1\), which yields the critical herdimmunity threshold^{8}
implying that a higher proportion needs to be vaccinated if the vaccine is not perfect^{8,9}. Since \(p^c\) must belong to the interval [0, 1], we set \(p^c\) equal to 1 if the value of (2) is larger than 1. We say a vaccine is effective if \(p^c<1\), and herd immunity is thus achievable.
Despite many advanced epidemiological models designed to explain the complex dynamics of infectious diseases^{10,11}, Eqs. (1) and (2) are the most applicable formulas in the history of vaccination. While (1) and (2) are based on an assumption of homogeneity of the infection in the population, new results show that (2) represents an upper bound for the herdimmunity threshold under a heterogeneity assumption^{9,12}. For example, it has been shown that in agestructured communities with heterogeneous activity groups, the diseaseinduced herdimmunity threshold would be lower than the one under the homogeneity assumption^{9}.
Achieving herd immunity is a multidimensional problem, depending not only on an individual’s decision but also the economic, environmental and societal conditions of the population. Considering all effective factors and overparameterising the epidemiological models not only makes them more complex but also makes it hard to extract plausible advice. To prevent overparameterization, we assume that \(\text {R}_0\) and E are constant and that only a single variant of the virus exists in the population, except in the penultimate section, which discusses a disease with multiple variants. Furthermore, in order to estimate the herdimmunity threshold, we assume that only pharmaceutical intervention is effective and that mass vaccination is the optimal way to control the disease. Henceforth, in this article, achieving herd immunity is equivalent to \(\text {R}_V\le 1\).
At the time of writing, multiple vaccines have received approval from governments, and their safety, efficacy, effectiveness and limitations have been analysed^{13,14,15,16}. There have been remarkable developments in RNAbased technologies to produce new vaccines against infectious diseases. It has been shown that the mRNA vaccines provide a longlasting immune response and they are potentially much safer than other vaccines. Such vaccines have the potential to be quickly manufactured and to become powerful tools against rapid pandemics. This new technology provides an opportunity to have more than one mRNA vaccine for any vaccinepreventable disease in the future^{17,18}. Therefore, extending Eqs. (1) and (2) is essential in order to better understand herd immunity under complex strategies and create a more reliable estimate of the herdimmunity threshold.
Suppose a policymaker is confronted with a pandemic and that multiple vaccines with different effectivenesses \(E_1, E_2, \ldots , E_k\) are available to control the spread. Using a single vaccine to vaccinate the entire population may not be the best strategy, because the available vaccine stockpile may be insufficient, and producing additional vaccines takes considerable time. In that situation, the policymaker must allocate different proportions of available vaccines to control the pandemic. However, if we want to use more than one vaccine, each with different effectiveness, for different proportions of the population, then Eqs. (1) and (2) are not valid. Here, we estimate the critical herdimmunity threshold for multivaccine allocation strategies. We also propose methods to allocate limited proportions of available vaccines in order to achieve herd immunity.
This article is organised as follows. First, we calculate the critical herdimmunity threshold under multivaccine strategies. Next, we discuss different methods to estimate optimal vaccineallocation strategies. We illustrate our results first for a multivaccine strategy for a single strain and secondly for the situation in Canada, with multiple variants and multiple vaccines. Finally, we conclude with a discussion.
Herd immunity and multiple vaccines
Assume that there are \(k>1\) effective vaccines for a single disease and that \(E_j\) is the effectiveness of the jth vaccine. Suppose that vaccinees have been selected randomly with the same probability among susceptible individuals. Policymakers can have many strategies to allocate vaccines based on various criteria. Although the most favourable vaccination strategy is using only a single vaccine with the largest effectiveness for the entire population, it is often impossible to receive enough of that vaccine in time, so a combination of different vaccines may be the optimal strategy. Therefore, policymakers have many possible choices to consider various proportions of each available vaccine to achieve herd immunity. We consider the least favourable vaccination strategy to be using only a single vaccine with the smallest effectiveness for the entire population.
Let \(p_j\) be the vaccineuptake proportion, a fraction of susceptible individuals who receive the jth vaccine, and let \(\text {S}\) be an allocation strategy with proportions \((p_1,p_2,\ldots , p_k)'\) such that \(\sum _{j=1}^kp_j\le 1\). Furthermore, assume that the susceptible populations receiving different vaccines are disjoint. The vaccination reproduction number with this multivaccine strategy is
which is the expected reproduction number in the population after vaccination. Our goal is to find a critical lower bound for the total proportion of the population, \(p_{_\text {S}}=\sum _{j=1}^k p_j\), that must be vaccinated under the vaccination strategy \(\text {S}\) to achieve \(\text {R}_V^{S}\le 1\). Obviously, the lower bound of \(p_{_\text {S}}\) is between critical herdimmunity thresholds for the least and most favourable strategies; i.e., it belongs to the interval
where \(p_j^c\) is the critical herdimmunity threshold if the jth vaccine were the only vaccine used. For any vaccination strategy, the critical threshold would be between the least or most favourable thresholds. The curves in Fig. 1 represent the herdimmunity threshold for different strategies as a function of the basic reproduction numbers when the effectiveness of each vaccines is constant. The black dashed curve and the black solid curve represent the mostfavourable and the leastfavourable vaccination strategies, respectively. Other strategies’ curves are between these two curves.
Consider a set of vaccines that are effective for small values of \(\text {R}_0\). Let \({\delta =\max _jp^c_j\min _jp_j^c}\) be the length of the interval (4), which represents the difference between the least favourable and the most favourable vaccination strategies. Figure 2 illustrates \(\delta\) as a function of the basic reproduction number, \(\text {R}_0\). Obviously, \(\delta\) is increasing for small values of \(\text {R}_0\). It is maximised at a point that the least favourable strategy is not effective anymore. The peak of the curve is at the last point of \(\text {R}_0\) where all vaccines are effective; i.e., the peak is at \(\text {R}_0=1/(1\min _{j=1}^k E_j)\). After this point, \(\delta\) decreases until the most favourable strategy is effective. Then, \(\delta\) is minimised to zero at \(\text {R}_0=1/(1\max _{j=1}^k E_j)\), which is the last point that the most favourable strategy is effective. After this point, \(\delta\) is constant and equal to zero; i.e., achieving vaccineinduced herd immunity is impossible under any combination of vaccines.
Suppose the policymaker has to select one of the available sets of vaccines for a single infectious disease. Each set defines a different scenario for vaccination and is illustrated by a different curve in Fig. 2. \(\delta\) can be interpreted as the effect of changing the share of a vaccine on the herdimmunity threshold. Usually, acceptance rates of vaccines are different, and their shares may change during a vaccination program, which results in a timedependent and variable herdimmunity threshold. A large value of \(\delta\) at the peak means that changing one vaccine’s share may substantially change the herdimmunity threshold. An important point is that the policymakers can propose effective multivaccine strategies if the reproduction number is close the peak of \(\delta\). On the other hand, improving a vaccination strategy is difficult for values of \(\text {R}_0\) when \(\delta\) is very small. Hence, replacing a vaccine with one that has a better effectiveness may not have a substantial effect on the herdimmunity threshold. Overall, we can say that the best set of vaccines is the one that takes its peak at a larger value of \(\text {R}_0\) and where \(\delta\) is small at the peak. An application of this plot is to compare vaccination programs in different countries. Each scenario in Fig. 2 can represent available vaccines for a given country.
Mean effectiveness of a multivaccine strategy
Estimating the smallest proportion of the population that must be vaccinated to achieve \(\text {R}_V^{S}\le 1\) is a challenging problem in controlling infectious diseases by vaccination. If we rewrite (3) as \(\text {R}_V^{\text {S}}=\text {R}_0\,(1p_S\bar{E}_S)\), where \(p_{_\text {S}}=\sum _{j=1}^k p_j\) is the total proportion of vaccinees, then \(\bar{E}_S\) is given by
which is a weighted mean of the effectiveness of all vaccines under the vaccination strategy \(\text {S}\). \(\bar{E}_{_\text {S}}\) measures the mean effectiveness of the vaccination strategy, which is different from the definition of the vaccine effectiveness. When we are dealing with multivaccine strategies, it shows how much a combination of vaccines would be efficient. Strategies with greater mean effectiveness are more promising candidates for an efficient impact.
For a singlevaccine strategy, \(\bar{E}_{_\text {S}}\) reduces to the effectiveness of a vaccine that is used for whole population. It is maximised for the most favourable and minimised for the least favourable vaccination strategy. When several vaccines are available, considering the vaccination of the population as a singlevaccine strategy whose effectiveness is the sample mean of all effectivenesses is equivalent to rewriting (3) as \(\text {R}_V^{\text {S}}=\text {R}_0\,(1p_S\bar{E})\), where \(\sum _{j=1}^kE_j/k\). This in valid only if the same proportion of all vaccines are available, which would not be a realistic assumption during most pandemics.
Relative effectiveness of a multivaccine strategy
In many cases, there are uncertainties in \(\text {R}_0\), and having a precise approximation is not possible. Assume that we can represent the uncertainties with an interval, \(\text {R}_0\in (a, b)\) such that a and b are known values. We define a measure of relative effectiveness for a strategy as \(\varepsilon ({\mathrm{S}})=1 {\displaystyle {\int }_a^b \mid p^c_{\mathrm {S}}p^c_{\mathrm {S}_m}\mid d\mathrm {R}_0}/{\displaystyle {\int }_a^b \mid p^c_{\mathrm {S}_\ell }p^c_{\mathrm {S}_m}\mid d\mathrm {R}_0}\), where \(\mathrm {S}_\ell\) and \(\mathrm {S}_m\) are the least and most favourable strategies, respectively. If all vaccines are effective for all values of \(\mathrm {R}_0\in (a, b)\), then it reduces to \(\varepsilon ({\mathrm{S}}) =1\frac{\mid 1/\bar{E}_{{\mathrm {S}}}1/\bar{E}_{{\mathrm{S}}_m}\mid }{\mid 1/\bar{E}_{{\mathrm {S}}_\ell }1/\bar{E}_{{\mathrm{S}}_m}\mid }\). Notice that \(\varepsilon (\mathrm {S})\in [0,1]\) and takes its maximum or minimum if the strategy is the most or least favourable strategy, respectively.
Herdimmunity threshold for a multivaccine strategy
Assume that \(p_j=\alpha _j\,p_j^c\), where \(\alpha _j\ge 0\) is a coefficient and \(p_j^c\) is the critical herdimmunity threshold for the jth vaccine if it were used for the whole population. Since the proportions \(p_j^c\) are known, given vaccine effectivenesses and the reproduction number of the disease, the vaccination strategy can be determined by \(\varvec{\alpha }=(\alpha _1,\ldots ,\alpha _k)'\). (Note that the coefficient of a vaccine is not the fraction of the population who receive that vaccine.) Hereafter, a vaccine’s coefficient is a key concept in the terminology of our discussion. The total proportion of the population that must be vaccinated under the strategy \(\mathrm {S}\) can be reformulated as
where (6) is a generalization of (2) and includes (2) as a special case if only a single vaccine exists. Furthermore, vaccineinduced herd immunity may be achieved if
or equivalently, if \(\mathrm {R}_V\le 1\).
The interpretation of the coefficients may be summarised as follows. In the absence of deriving a closed formula for the herdimmunity threshold from the equation \(R_v=R_0(1\sum _j p_j E_j)\), we formulate the problem using \(p_j=\alpha _j p^c_j=\alpha _j (11/R_0) /E_j\), or equivalently \(\alpha _j=p_j/p_j^c\) (where \(p_j^c\) is a function of known values \(R_0\) and \(E_j\)). This formulation leads to the condition \(\sum _j\alpha _j\ge 1\) for achieving the herd immunity. Hence \(p_j\) is the weight of the jth vaccine with corresponding coefficient \(\alpha _j\).
When a vaccination plan changes
When a vaccination process starts, it can be stopped at any time because of new evidence on the vaccine’s side effects or any other reasons. Many social, political and economical factors may have substantial effects on vaccination progress. Assume that a vaccination process started with k vaccines under a strategy \(\mathrm {S}_{(1)}\) with proportions \((p_{11}^*,\ldots ,p_{1k}^*)'\), but after a while it has stopped. After or without a pause, the vaccination process will continue with adding, removing or replacing some vaccines. Assume that the second step of the vaccination contains \(k'\) vaccines with effectivenesses \(E_{21},\ldots , E_{2k'}\) under a new strategy \(\mathrm {S}_{(2)}\) with proportions \((p_{21},\ldots , p_{2k'})'\). The vaccination reproduction number is \(\mathrm {R}_V^{\mathrm {S}_{(2)}} =\mathrm {R}_0\,\big (1q\sum _{j=1}^{k'}p_{2j}E_{2j}\big )\), where \(q=\sum _{j=1}^kp_{1j}^*E_{1j}\) is an approximation of the proportion of individuals who are immunised by vaccination in the first step. Note that q can be considered as a proportion of individuals who are immunised naturally by infection if the disease provides longterm immunity and if people who are infected are excluded from the vaccination program. Setting \(p_{2j}=\alpha _{2j}\,p^c_{2j}\) implies that herd immunity can be reached if
which leads to a smaller level of vaccination coverage in order to achieve herd immunity. Note that this result can be extended to cases where the vaccination process has more than two steps.
Strategy estimation
When an infection exists in the population and vaccination is the main tool to stop spreading it, any failure to stockpile sufficient doses of the vaccine can have drastic consequences on public health. Hence, finding an optimal vaccination strategy that leads to herd immunity is an important step in fighting the disease. In this section, we introduce methods to estimate the unknown vaccine coefficients \(\alpha _j\) for \(j=1,\ldots , k\) such that k is a known integer and \(E_1, \ldots , E_k\) are known effectivenesses of available vaccines.
Rankingbased method
Assigning coefficients to vaccines, calculated with respect to some objective criteria, is a simple way to estimate the optimal vaccination strategies. Let u be a positive variable that represents the utility of the vaccines, such that if \(u_i<u_j\), then the jth vaccine must have a larger coefficient than the ith vaccine in a vaccination strategy. With this ordering, we can define
where n is a shrinkage parameter that takes nonnegative values. By controlling n, some coefficients can be shrunk to zero, and their corresponding vaccines are removed from the vaccination program if desired. Here, we introduce some rankingbased methods:

(A)
When there is no prior information to sort available vaccines, a simple strategy with equal coefficients for all vaccines can be used; i.e., \(\alpha _j= {1}/{k}\). This strategy can be used when there is no meaningful difference between vaccines with respect to decision components. This strategy is equivalent to the case where we simply apply a classic singlevaccine strategy whose effectiveness is the average of the effectivenesses of all the available vaccines.

(B)
When effectiveness is the most important criteria to sort vaccines, we may set \(u_j=E_j\). This strategy may be used for the case that all vaccines are effective but variation between their effectivenesses is large.

(C)
If reducing side effects of the vaccines is the most important criteria, we must maximise the proportion of individuals who remain both unvaccinated and uninfected. This proportion is known as herd effect^{19}. In this situation, we can set \(u_j=1p_j^c\) where \(p_j^c\) is the herdimmunity threshold for the jth vaccine. Under this strategy, if a vaccine is not effective and we cannot stop the spread of the disease by using only that specific vaccine, its coefficient will be zero.

(D)
A strategy can be based on some criteria such as the length of vaccineinduced immunity, vaccine availability, vaccination cost, etc. If we assume that vaccines are ranked based on multiple criteria where \(r_{j}\) is the rank of the jth vaccine such that larger ranks are more important than smaller ones, then we can set \(u_j=r_j\). This strategy would be the best choice when many factors must be considered to order the vaccines.
Optimising herd effect and vaccination cost
Usually vaccinating the entire population is impossible, since immunocompromised people, pregnant women or young children often cannot be vaccinated. Therefore, the ideal vaccination policy would be the one with the lowest herdimmunity threshold, which allows the most people to avoid the infection without vaccination. Optimisation techniques can be used to estimate the optimal strategy with the appropriate proportion of vaccines, maximising the benefits while minimising the costs. Optimal vaccination strategies to fight against specific disease for a single vaccine have been discussed^{20,21}. Let \(\mathrm {S}\) be a vaccination strategy with unknown coefficients \(\varvec{\alpha }=(\alpha _1,\ldots , \alpha _k)'\) such that \(\sum _{j=1}^k\alpha _j\ge 1\). A reasonable criterion for determining the optimal vaccineallocation strategy is achieving herd immunity such that the number of individuals who escape the infection without vaccination is maximised, which indirectly decreases the side effects of the vaccine. This is equivalent to minimising \(\sum _{j=1}^k {\alpha _j}/{E_j}\), such that \(\sum _{j=1}^k\alpha _j\ge 1\).
Optimising the total cost of vaccination is an important consideration for policymakers. Vaccination can have a major effect on developing countries, but health budgets are often limited. Therefore, cost effectiveness is a necessary part of mathematical modelling of the vaccination process. Let N be the size of the population of interest, \(c_j\) be the cost of the jth vaccine per individual, for \(j=1, \ldots , k\). Then a lower bound for the total cost of vaccination under the vaccination policy \(\mathrm {S}\) is \(N(1\frac{1}{\mathrm {R}_0})\sum _{j=1}^k\alpha _j {c_j}/{E_j}\). Assuming N and \(\mathrm {R}_0\) are constant, one can estimate the vaccine weights by minimising the cost and maximizing the herd effect simultaneously, which yields the convex optimisation problem
where all \(c_j\) and \(E_j\) are known constant values.
Example 1
Consider a situation of an infectious disease spreading in the population. Here, we explain how our proposed method can be used to allocate specific proportion of available vaccines to reach the herd immunity.
Assume that \(\mathrm {R}_0=2.5\) and six vaccines \(V_1,\ldots , V_6\) are available with different efficiencies \(E_1=0.80\), \(E_2=0.85\), \(E_3=0.90\), \(E_4=0.94\), \(E_5=0.95\), and \(E_6=0.95\); also, assume \(\mathrm {S}_{V_j}\) represents a singlevaccine strategy such that only the jth vaccine is used for the entire population; \(\mathrm {S}_{A}^n\), \(\mathrm {S}_{B}^n\), \(\mathrm {S}_{C}^n\), and \(\mathrm {S}_{D}^n\) are multivaccine strategies based on a rankingbased method, such that n can take values 5, 10 or 20. For the rankbased method, vaccines are ranked based on ascending order of their effectivenesses such that \(r_i=i\) and the average of orders is used when ties exist and some effectivenesses are equal. All the information is presented in Table 1. The last three columns of the table represent \(p^c_{\mathrm {S}}\), the total proportion of individuals who must be vaccinated under each strategy, relative effectiveness of strategies, \(\varepsilon (\mathrm {S})\), and mean effectiveness of strategies, \(\bar{E}_{\mathrm {S}}\), respectively. The table body elements are \(\alpha _j\) for \(j=1,\ldots , 6\).
Table 1 is divided into two parts: the first represents singlevariance strategies, while the second describes multivaccine strategies. The first part consists of six singlevaccine strategies such that \(p_j^c=P_{\mathrm {S}_{V_j}}^c\) for \(j=1, \ldots , 6\). Therefore, \(p_1^c=0.75\), \(p_2^c=0.706\), \(p_3^c=0.667\), \(p_4^c=0.638\), \(p_5^c=0.632\) and \(p_6^c=0.632\) are critical vaccine thresholds for \(V_1,\ldots , V_6\), respectively. The least favourable strategy \(\mathrm {S}_{V_1}\) implies that at least \(75\%\) of the population must be vaccinated while the most favourable strategy needs only \(63.2\%\) of the population to be vaccinated.
Note that an infinite number of multivaccine strategies can be proposed. However, only a few strategies are discussed, based on estimation methods described in the rankingbased method. There are four categories of multivaccine strategies in the second part of Table 1:

The first category of the second part is \(\mathrm {S}_A\), where uniform coefficients are assigned to the vaccines. It leads to vaccination coverage equal to 0.67, which is equivalent to vaccinating the entire population with a single vaccine with effectiveness of \(E=0.90\). Since \(\alpha _j=1/6\approx 0.167\) and \(p_j=\alpha _jp^c_j\), then \(p_1\approx 0.125\), \(p_2\approx 0.118\), \(p_3\approx 0.111\), \(p_4\approx 0.107\), \(p_5\approx 0.106\) and \(p_6\approx 0.106\) are proportions of the population that will be vaccinated by \(V_1,\ldots , V_6\) to reach herdimmunity.

The second category is based on effectiveness values. By increasing the shrinkage parameter, the coefficients approaches zero slower than strategies in Categories C and D. For example, when \(n=20\), Strategy \(\mathrm {S}_B^{20}\) proposes to use all vaccines, while \(\mathrm {S}_C^{20}\) and \(\mathrm {S}_D^{20}\) remove some vaccines from the vaccination program.

Category C of multivaccine strategies is based on herd immunity and assigns the largest coefficient to a vaccine that lets more people escape the infection without vaccination. Strategies in Category C can be close to some strategy in Category B. For example, strategies \(\mathrm {S}_C^{10}\) and \(\mathrm {S}_B^{20}\) suggest similar coefficients for vaccines.

The last category is based on ranking vaccines. The shrinkage parameter reduces the coefficients faster than Categories B and C. For example, \(\mathrm {S}_D^{20}\) assigns coefficients only to vaccines with largest effectiveness.

For a constant value of the shrinkage parameter \(n=n_0\), \(\mathrm {S}_D^{n_0}\) is better than \(\mathrm {S}_C^{n_0}\) and \(\mathrm {S}_C^{n_0}\) is better than \(\mathrm {S}_B^{n_0}\), according to their mean and relative effectivenesses. When n is extremely large, all strategies in Categories B, C and D approach \(\mathrm {S}_D^{20}\), which is the most favourable strategy.

Note that relative effectiveness \(\varepsilon (\mathrm {S})\) represents the power of a possible combination of available vaccines with respect to the least and most favourable strategies; its small value does not imply weakness of the strategy to stop the infection. For example, the relative effectiveness is always small for strategies that put large weight on the least favourable vaccine regardless of its effectiveness. Hence, both mean and relative effectivenesses together show the power of a multivaccine strategy to contain the infection.
Making decisions for vaccine allocation is a complex problem and depends on multifactor priorities. Table 1 proposes plausible strategies that satisfy logistic, financial and social conditions of the publichealth system. Overall, all strategies in Categories B, C, and D are better than \(\mathrm {S}_A\), which is equivalent to a singlevaccine strategy whose effectiveness is the sample mean of all effectivenesses.
When several variants exist
When several variants of a virus exist in the population, some adjustments are needed to estimate the herdimmunity threshold. Assume that only a single vaccine exists to control spreading a virus with m variants and that \(\pi _1,\ldots , \pi _m\) are the proportions of each variant, such that \(\sum _{j=1}^m \pi _j=1\). In addition, \(E_1, \ldots , E_m\) are effectivenesses of the vaccine against variants. Under this setting, the vaccination reproduction number is \(\bar{{\mathrm{R}}}_V=\bar{{\mathrm{R}}}_0\,(1p\bar{E})\), where \(\bar{\mathrm {R}}_0=\sum_{j=1}^m\pi_j\mathrm{R}_0^j\) is the average reproduction number for all variants and \(\bar{E}=\sum _{j=1}^m \pi _j E_j\) is the average effectiveness of the vaccine against the virus in the population. Note that estimating \(\bar{E}\) would be a challenging problem if adequate and precise clinical tests to diagnose different variants are not done with respect to spatial distribution of the disease, which could lead to biased estimation of the prevalence of different variants. In addition, distribution of variants is timedependent; consequently, \(\bar{E}\) will change over time.
If k vaccines are available for a virus with m variants, we assume that \(E_{ij}\) is the effectiveness of the ith vaccine against the jth variant. Therefore, the vaccination reproduction number would be \(\bar{{\mathrm{R}}}_V=\bar{{\mathrm{R}}}_0\,(1\sum _{j=1}^k p_i \bar{E}_i)\), where \(\bar{E}_i=\sum _{j=1}^m \pi _j E_{ij}\) is the average effectiveness of the ith vaccine against the virus. Herd immunity may thus be achieved if
where \(p_i=\alpha_i[(11/R_0)/E_{i1}]\) and \(E_{i1}\) is the effectiveness of the ith vaccine against the original strain of the virus. This implies that when new emerging variants reduce the effectiveness of the vaccines, more vaccines are needed to reach herd immunity.
Example 2
In this example, we show how to estimate the vaccination reproduction number when several variants with different prevalence fractions exist in the population and multiple vaccines are used to control the spread of the disease. Multiple studies have evaluated the effectiveness of different vaccines against important variants of COVID19^{22}. Realworld data on vaccine effectiveness against all variants recognized by the World Health Organization (WHO) are currently limited. New emerging variants of concern (VOC)—currently including Alpha (B.1.1.7), Beta (B.1.351), Gamma (P.1) and Delta (B.1.617.2)—are more contagious and can potentially increase disease severity and decrease vaccine effectiveness^{22,23,24,25,26}.
Vaccination against symptomatic infection of COVID19 in 2020–2021 has been deployed using BNT162b2 (PfizerBioNTech), mRNA1273 (Moderna) and ChAdOx1 (AstraZeneca) vaccines^{26}. The effectiveness of these vaccines against VOC has been estimated by Nasreen et al.^{26}. Estimates of the effectiveness of vaccines against infection with partial and full vaccination in Canada is reported in Table 2. However, due to insufficient data for Moderna and AstraZeneca vaccines, effectiveness of fully vaccinated cases was not estimated for these vaccines. Therefore, we estimated missing values and used other sources to provide estimates for data that we need. We make the following observations about Table 2.

The second column represents the type of vaccine. In some cases, the vaccine type was not reported or individuals used mixed vaccines; i.e., first and second doses were not from the same vaccine.
The third column presents cumulative percentages of people who had received a COVID19 vaccine by July 31, 2021. This information is available on website of the government of Canada at https://healthinfobase.canada.ca.
The fourth through seventh columns represent effectiveness of vaccines against the original strain and VOC. Since both Beta and Gamma variants share common mutations^{26} and insufficient numbers of them were sequenced, they are classified together.

A report from Public Health Ontario, https://www.publichealthontario.ca, provides a list of the reported effectiveness measures for vaccines around the world. In general, vaccine effectiveness for preventing symptomatic infection 3–4 weeks after receiving a single dose is between 60% and 80% for the PfizerBioNTech and Moderna, and between 60% and 70% for AstraZeneca. Since the effectiveness of vaccines decreases against VOC, we estimate the effectiveness of PfizerBioNTech, Moderna and AstraZeneca for a single dose against the original strain at 80%, 80% and 70%, respectively.

The Public Health England vaccine effectiveness report from March 2021, https://assets.publishing.service.gov.uk, shows that vaccine effectiveness against infection is about 85%, seven or more days after the second dose.

Lopez Bernal et al. showed that vaccine effectiveness after two doses of the AstraZeneca vaccine is 74% against the alpha variant and 67% against the delta variant^{25}. We used this information to estimate the missing values.

The final column represents average effectiveness of each vaccine against all listed variants, which is calculated as a weighted mean of effectiveness measures; i.e., \(\bar{E}_i=\sum _{j}\pi _jE_{ij}\). The weights are prevalence of variants, which is represented in the bottom row of Table 2. This row presents cumulative prevalence of VOC as listed in the Canada mutation report, dated August 15, 2021. It is available at https://outbreak.info, which provides results of screening tests for mutations and wholegenome sequencing to assign COVID19 lineage. The cumulative prevalence is the ratio of the genetic sequences containing the lineage or mutations to all genetic sequences collected since the identification of lineage or mutations in that location.

The fourth row of singledose and twodose vaccination show “Vaccines not reported”, which represent a group of people who are vaccinated by PfizerBioNTech, Moderna or AstraZeneca, but the name of their vaccines are missing from the reported data. These account for 3.33% and 13.2% of people take singledose and twodose vaccines, respectively. We have estimated their corresponding rows by a weighted mean of effectiveness of all the vaccines with respect to their cumulative percentages in third column.

The vaccine effectiveness of Moderna for full vaccination was estimated only against Alpha but was missing for other variants. We used the effectiveness of Moderna after first dose to estimate the missing values.

Note that only 0.1% and 0.02% of people were partially and fully vaccinated using other vaccines. Therefore, we ignored this category.

In order not to underestimate the reproduction number of the virus, we estimate missing values with smallest similar data. For example, in full vaccination, effectiveness of AstraZeneca was available only for the original variant, so we estimate its effectiveness against VOC with the same value.

Mixed vaccination was approved by the government of Canada, and a proportion of people used different vaccines for their first and second doses. Since it was not clear which vaccines were mixed together, we assumed that most of them used AstraZeneca for the first dose and Moderna or PfizerBioNTech for the second dose. Therefore, to prevent overestimating the missing effectiveness values, we estimated the effectiveness of the mixed vaccines with the smallest effectiveness of those three vaccines for full vaccination.

Empirical research shows substantial increases in estimated reproduction numbers of VOC, compared to the original strain. These increases are 29%, 25%, 38% and 97%, respectively, for Alpha, Beta, Gamma and Delta variants. We considered an average of 30% increase for Beta/Gamma^{27}.
We can thus estimate the vaccination reproduction number using our proposed method. If we set \(\mathrm {R}_0=2.5\) for the original strain, we have \({\mathrm{R}}_V = {\mathrm{R}}_0\, ( 1\sum _{_{Partial}} p_i E_{i1} \sum _{_{Full}} p_i E_{i1} )=0.8817\) when only the original variant exists in the population. However, if several variants coexist and we modify the reproduction number with respect to their prevalence distribution and effectiveness of vaccines against them according to Table 2, then \(\bar{{\mathrm{R}}}_V=\bar{{\mathrm{R}}}_0\, ( 1\sum _{_{Partial}} p_i \bar{E}_{i} \sum _{_{Full}} p_i \bar{E}_{i} )=1.17\) is our modified reproduction number. It follows that \(\bar{\mathrm {R}}_V/ {\mathrm {R}}_V\approx 1.32\); hence if we do not modify our calculations in presence of new emerging variants that are potentially resistant against vaccines, we would underestimate the reproduction number, which may lead us to biased decisions. Note that \(\bar{\mathrm {R}}_V/ \mathrm {R}_V\) may change substantially if the prevalence of new emerging vaccineresistant variants increases. This example illustrates how multiple variants may affect the vaccination reproduction number.
Discussion
We have analysed the complexity of herd immunity in the presence of multiple vaccines and multiple variants and proposed an update to the herdimmunity threshold when multiple vaccines are used or when several variants of the virus coexist. This update considers not only distribution of the vaccines but also the prevalence of variants and effectiveness of vaccines against them. Our contributions bring several interesting insights to the literature of herd immunity and vaccineallocation strategy. First, we calculated the herdimmunity threshold for combinations of multiple vaccines and showed that it depends on not only the reproduction number and vaccine effectiveness but also the proportion of vaccinees for all vaccines. Moreover, we determined which combination of vaccines would lead to herd immunity. Second, we showed that when several vaccines are available, simplifying the problem and considering the vaccination of the population as a singlevaccine strategy whose effectiveness is the sample mean of all effectivenesses, \(\sum _{j=1}^kE_j/k\), is not ideal, because many multivaccine strategies with smaller herdimmunity thresholds can be proposed. Third, our proposed modification of the herdimmunity threshold enables us to search for optimal combination of vaccines in order to reach herd immunity. We described methods to search for optimal vaccineallocation strategies based on different priorities. Fourth, we showed that vaccineallocation policies are equivalent to simple convex optimisation problems if we wish to maximise the herd effect and/or minimise the cost of a vaccination strategy. Fifth, our definition of the mean effectiveness, which is a novel measure for describing efficiency of a vaccination program, can be used to compare different vaccination strategies. Finally, we extended the main results of this work to the situation where multiple variants of concern exist in the population. We determined conditions for reaching herd immunity in the presence of new emerging variants of concern. We have shown that several coexisting variants may change the vaccination effectiveness and have illustrated it with a Canadian example. It follows that the herdimmunity threshold must be updated with respect to the prevalence distribution of variants and effectiveness of vaccines against them.
Vaccine hesitancy is one of the most important problems in achieving herd immunity. The WHO has listed vaccine hesitancy as one of ten threats to global health. Vaccine hesitancy is defined as a delay in acceptance or refusal of vaccines despite the availability of vaccination services^{28}. Recent investigations report different vaccineacceptance rates for available vaccines for COVID19^{29,30}. Therefore, vaccineacceptance rates or vaccinehesitancy rates are potential components that can be used to optimise the vaccine’s coefficients in order to introduce an optimal vaccination strategy that minimises overall vaccine hesitancy.
Having more vaccines with large efficacies increases the strength of any vaccination program. However, the speed of vaccination should not be underrated in pandemic management. Containing the disease in the shortest possible time is an ideal goal, which needs to use the maximum amount of available vaccines. Therefore, having more available vaccines would be considered a component of an optimal strategy.
If we assume that an individual can be immunised after infection for a short period of time, then the vaccination reproduction number is timedependent and can be modified as \({\mathrm{R}}_V(t) = {\mathrm{R}}_0(1q(t)\sum_{j=1}^kp_jE_j)\), where \(q(t)\) is the proportion of the population that is still immunised naturally at time t, which can be estimated empirically from antibody test results or from related epidemiological models. This implies that, with natural immunisation after infection, herd immunity may be achieved faster than predicted by a vaccineinduced herdimmunity threshold. In contrast, we have shown that when multiple variants of the virus are emerging in the population, more vaccines are needed to reach the herdimmunity level. Therefore, a generalised condition for reaching herd immunity is \(\sum_{j=1}^k\alpha_j(\bar E_i/E_{i1})\geq1q(t)/(11/\mathrm{R}_0)\), which includes (8) and (11) as special cases. Hence it is reasonable to expect to reach herd immunity at larger vaccination levels when the probability of reinfection is significantly large.
Our results have some limitations, which should be acknowledged. The first is that there is no measurement error in the estimates of R_{0} or E. the second is that they are constant across all communities in the population, which implies that the population is homogeneous; this is not true in agestructured populations, for instance. The third is that they are not timedependent. Further work is required to generalise our results in these circumstances.
Vaccination at a level larger than the herdimmunity threshold does not imply that the infection can be stopped if vaccination has a different pattern from the infection. For example, a vaccination policy based on age prioritisation can save more lives if age is the only heterogeneous characteristic of the population and infection is uniformly distributed spatially. Without considering the spatial distribution of the infection, the vaccination process can waste valuable time and cause the disease to spread faster in the population. In an effective vaccination policy, spatial hotspots must take priority, and vaccines must therefore be distributed with the same spatial pattern as the infection itself. Only if proactive measures are taken will we be able to prevent the most infections and save the most lives using all possible vaccines at our disposal.
References
 1.
Haffajee, R. L. & Mello, M. M. Thinking globally, acting locally—The US response to COVID19. N. Engl. J. Med. 382, e75 (2020).
 2.
Lu, N., Cheng, K.W., Qamar, N., Huang, K.C. & Johnson, J. A. Weathering COVID19 storm: Successful control measures of five Asian countries. Am. J. Infect. Control 48, 851–852 (2020).
 3.
Jeyanathan, M. et al. Immunological considerations for COVID19 vaccine strategies. Nat. Rev. Immunol. 20, 615–632 (2020).
 4.
Metcalf, C. J. E., Ferrari, M., Graham, A. L. & Grenfell, B. T. Understanding herd immunity. Trends Immunol. 36, 753–755 (2015).
 5.
Anderson, R. M. & May, R. M. Vaccination and herd immunity to infectious diseases. Nature 318, 323–329 (1985).
 6.
Delamater, P. L., Street, E. J., Leslie, T. F., Yang, Y. T. & Jacobsen, K. H. Complexity of the basic reproduction number (R_{0}). Emerg. Infect. Dis. 25, 1 (2019).
 7.
Anderson, R. M. & May, R. M. Immunisation and herd immunity. Lancet 335, 641–645 (1990).
 8.
Fine, P., Eames, K. & Heymann, D. L. Herd immunity: A rough guide. Clin. Infect. Dis. 52, 911–916 (2011).
 9.
Britton, T., Ball, F. & Trapman, P. A mathematical model reveals the influence of population heterogeneity on herd immunity to SARSCoV2. Science 369, 846–849 (2020).
 10.
Hethcote, H. W. The mathematics of infectious diseases. SIAM Rev. 42, 599–653 (2000).
 11.
Wang, Z. et al. Statistical physics of vaccination. Phys. Rep. 664, 1–113 (2016).
 12.
Gomes, M. G. M. et al. Individual variation in susceptibility or exposure to SARSCoV2 lowers the herd immunity threshold. Preprint at medRxiv https://doi.org/10.1101/2020.04.27.20081893 (2020).
 13.
Voysey, M. et al. Safety and efficacy of the ChAdOx1 nCoV19 vaccine (AZD1222) against SARSCoV2: An interim analysis of four randomised controlled trials in Brazil, South Africa, and the UK. Lancet (2020).
 14.
Knoll, M. D. & Wonodi, C. OxfordAstraZeneca COVID19 vaccine efficacy. Lancet 397, 72–74 (2021).
 15.
Burki, T. K. The Russian vaccine for COVID19. Lancet Respir. Med. 8, e85–e86 (2020).
 16.
Polack, F. P. et al. Safety and efficacy of the bnt162b2 mRNA COVID19 vaccine. N. Engl. J. Med. 383, 2603–2615 (2020).
 17.
Pardi, N., Hogan, M. J., Porter, F. W. & Weissman, D. mRNA vaccines—A new era in vaccinology. Nat. Rev. Drug Discov. 17, 261 (2018).
 18.
Mascola, J. R. & Fauci, A. S. Novel vaccine technologies for the 21st century. Nat. Rev. Immunol. 20, 87–88 (2020).
 19.
John, T. J. & Samuel, R. Herd immunity and herd effect: New insights and definitions. Eur. J. Epidemiol. 16, 601–606 (2000).
 20.
Lemmens, S., Decouttere, C., Vandaele, N. & Bernuzzi, M. A review of integrated supply chain network design models: Key issues for vaccine supply chains. Chem. Eng. Res. Des. 109, 366–384 (2016).
 21.
Duijzer, L., Van Jaarsveld, W., Wallinga, J. & Dekker, R. Doseoptimal vaccine allocation over multiple populations (tech. rep.). Econometric Institute, Erasmus School of Economics (2015).
 22.
Imai, N. et al. Interpreting estimates of coronavirus disease 2019 (covid19) vaccine efficacy and effectiveness to inform simulation studies of vaccine impact: A systematic review. Wellcome Open Res. 6, 185 (2021).
 23.
Sheikh, A., McMenamin, J., Taylor, B. & Robertson, C. Sarscov2 delta voc in scotland: Demographics, risk of hospital admission, and vaccine effectiveness. Lancet 397, 2461–2462 (2021).
 24.
AbuRaddad, L. J., Chemaitelly, H. & Butt, A. A. Effectiveness of the bnt162b2 covid19 vaccine against the b. 1.1. 7 and b. 1.351 variants. N. Engl. J. Med. 385, 187–189 (2021).
 25.
Lopez Bernal, J. et al. Effectiveness of covid19 vaccines against the B.1.617. 2 (Delta) variant. N. Engl. J. Med. 385, 585–594 (2021).
 26.
Nasreen, S. et al. Effectiveness of mRNA and ChAdOx1 COVID19 vaccines against symptomatic SARSCoV2 infection and severe outcomes with variants of concern in Ontario. Preprint at medRxiv https://doi.org/10.1101/2021.06.28.21259420 (2021).
 27.
Campbell, F. et al. Increased transmissibility and global spead of SARSCOV1 variants of concern as at June 2021. Eurosurveillance 26, 2100509 (2021).
 28.
MacDonald, N. E. et al. Vaccine hesitancy: Definition, scope and determinants. Vaccine 33, 4161–4164 (2015).
 29.
Wouters, O. J. et al. Challenges in ensuring global access to COVID19 vaccines: Production, affordability, allocation, and deployment. Lancet 397, 1023–1034 (2021).
 30.
Schwarzinger, M., Watson, V., Arwidson, P., Alla, F. & Luchini, S. Covid19 vaccine hesitancy in a representative workingage population in France: A survey experiment based on vaccine characteristics. Lancet Public Health 6(4), e210–e211 (2021).
Acknowledgements
The authors are grateful to two anonymous reviewers whose comments greatly improved the manuscript. S.R.S? is supported by NSERC Alliance and Discovery Grants.
Author information
Affiliations
Contributions
I.Y. conceptualized the problem, codesigned the study and cowrote the manuscript, O.M. codesigned the study, S.R.S.? cowrote and edited the manuscript. All authors reviewed the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher's note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Yadegari, I., Omidi, M. & Smith, S.R. The herdimmunity threshold must be updated for multivaccine strategies and multiple variants. Sci Rep 11, 22970 (2021). https://doi.org/10.1038/s41598021000832
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598021000832
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.