The herd-immunity threshold must be updated for multi-vaccine strategies and multiple variants

Several vaccines with different efficacies and effectivenesses are currently being distributed across the world to control the COVID-19 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 vaccine-induced herd immunity by vaccinating a fraction of the population. The classic vaccine-induced herd-immunity 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 multi-vaccine strategies for a disease or multiple variants, potentially leading policymakers to suboptimal vaccine-allocation 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 single-vaccine strategy whose effectiveness is the sample mean of all effectivenesses would not be ideal, because many multi-vaccine strategies with a smaller herd-immunity threshold can be proposed. We show that the herd-immunity threshold may vary due to changes in vaccine-uptake 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 herd-immunity threshold must be updated not only when multi-vaccine strategies are used but also when multiple variants coexist in the population.

Following the rapid worldwide spread of the coronavirus SARS-CoV-2, 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 public-health 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 herd-immunity 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, 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 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 R 0 is less than one and to increase if 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

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 j-th 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 vaccine-uptake proportion, a fraction of susceptible individuals who receive the j-th vaccine, and let S be an allocation strategy with proportions (p 1 , p 2 , . . . , p k ) ′ such that k j=1 p j ≤ 1 . Furthermore, assume that the susceptible populations receiving different vaccines are disjoint. The vaccination reproduction number with this multi-vaccine strategy is www.nature.com/scientificreports/ 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 S = k j=1 p j , that must be vaccinated under the vaccination strategy S to achieve R S V ≤ 1 . Obviously, the lower bound of p S is between critical herd-immunity thresholds for the least and most favourable strategies; i.e., it belongs to the interval where p c j is the critical herd-immunity threshold if the j-th 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 herd-immunity 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 most-favourable and the least-favourable vaccination strategies, respectively. Other strategies' curves are between these two curves.
Consider a set of vaccines that are effective for small values of R 0 . Let δ = max j p c j − min j p c j be the length of the interval (4), which represents the difference between the least favourable and the most favourable vaccination strategies. Figure 2 illustrates δ as a function of the basic reproduction number, R 0 . Obviously, δ is increasing for small values of 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 R 0 where all vaccines are effective; i.e., the peak is at R 0 = 1/(1 − min k j=1 E j ) . After this point, δ decreases until the most favourable strategy is effective. Then, δ is minimised to zero at , which is the last point that the most favourable strategy is effective. After this point, δ is constant and equal to zero; i.e., achieving vaccine-induced 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. δ can be interpreted as the effect of changing the share of a vaccine on the herd-immunity threshold. Usually, acceptance rates of vaccines are different, and their shares may change during a vaccination program, which results in a time-dependent and variable herd-immunity threshold. A large value of δ at the peak means that changing one vaccine's share may substantially change the herd-immunity threshold. An important point is that the policymakers can propose effective multi-vaccine strategies if the reproduction number is close the peak of δ . On the other hand, improving a vaccination strategy is difficult for values of R 0 when δ is very small. Hence, replacing a vaccine with one that has a better effectiveness may not have a substantial effect on the herd-immunity threshold. Overall, we can say that the best set of vaccines is the one that takes its peak at a larger value of R 0 and where δ 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 multi-vaccine strategy. Estimating the smallest proportion of the population that must be vaccinated to achieve R S V ≤ 1 is a challenging problem in controlling infectious diseases by vac-(4) min j p c j , max j p c j , Figure 1. Illustration of the herd-immunity thresholds for four different vaccine-allocation strategies as a function of R 0 when k = 6 vaccines are available and effectivenesses belong to the interval (0.7, 0.975): the black dashed curve represents the most favourable strategy, and the black solid curve is the least favourable strategy; all vaccine-allocation strategies will be between these two curves; the red dotted curve and the blue dashed-dotted curve are two different examples of vaccines strategies, represented by S A and S B , respectively. The strategy S A is equivalent to a single-vaccination strategy whose effectiveness is the average of effectivenesses. We can also propose strategies like S B that are closer to the most favourable strategy.
where p S = k j=1 p j is the total proportion of vaccinees, then Ē S is given by which is a weighted mean of the effectiveness of all vaccines under the vaccination strategy S . Ē S measures the mean effectiveness of the vaccination strategy, which is different from the definition of the vaccine effectiveness. When we are dealing with multi-vaccine 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 single-vaccine strategy, Ē 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 single-vaccine strategy whose effectiveness is the sample mean of all effectivenesses is equivalent to rewriting (3) as R S V = R 0 (1 − p SĒ ) , where k j=1 E 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 multi-vaccine strategy. In many cases, there are uncertainties in R 0 , and having a precise approximation is not possible. Assume that we can represent the uncertainties with an interval, R 0 ∈ (a, b) such that a and b are known values. We define a measure of relative effectiveness for a strategy as where S ℓ and S m are the least and most favourable strategies, respectively. If all vaccines are effective for all values of R 0 ∈ (a, b) , then it reduces to ε(S) Notice that ε(S) ∈ [0, 1] and takes its maximum or minimum if the strategy is the most or least favourable strategy, respectively.
Herd-immunity threshold for a multi-vaccine strategy. Assume that p j = α j p c j , where α j ≥ 0 is a coefficient and p c j is the critical herd-immunity threshold for the j-th vaccine if it were used for the whole population. Since the proportions p c j are known, given vaccine effectivenesses and the reproduction number of the disease, the vaccination strategy can be determined by α = (α 1 , . . . , α 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 S can be reformulated as Figure 2. Illustration of the length of the interval (4) as a function of R 0 for three different scenarios: the blue solid curve represents vaccine effectivenesses between 0.7 and 0.95; the red dashed curve represents effectivenesses between 0.7 and 0.9; the black dotted curve represents effectivenesses between 0.85 and 0.96. The first and second scenarios have a peak at the same point with different heights. However, the third scenario has a peak at a larger R 0 value with a smaller height at its peak. Hence the third set of vaccines is effective for a larger range of reproduction numbers, and the difference between the least and the most favourable strategies is smaller than other scenarios.  (6) is a generalization of (2) and includes (2) as a special case if only a single vaccine exists. Furthermore, vaccine-induced herd immunity may be achieved if The interpretation of the coefficients may be summarised as follows. In the absence of deriving a closed formula for the herd-immunity threshold from the equation where p c j is a function of known values R 0 and E j ). This formulation leads to the condition j α j ≥ 1 for achieving the herd immunity. Hence p j is the weight of the j-th vaccine with corresponding coefficient α 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 S (1) with proportions (p * 11 , . . . , 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 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 long-term immunity and if people who are infected are excluded from the vaccination program. Setting p 2j = α 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 α j for j = 1, . . . , k such that k is a known integer and E 1 , . . . , E k are known effectivenesses of available vaccines.
Ranking-based 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 j-th vaccine must have a larger coefficient than the i-th vaccine in a vaccination strategy. With this ordering, we can define where n is a shrinkage parameter that takes non-negative 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 ranking-based 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., α 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 single-vaccine 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 = 1 − p c j where p c j is the herd-immunity threshold for the j-th 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.
k j=1 α j ≥ 1, www.nature.com/scientificreports/ (D) A strategy can be based on some criteria such as the length of vaccine-induced 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 j-th 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 herd-immunity 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 S be a vaccination strategy with unknown coefficients α = (α 1 , . . . , α k ) ′ such that k j=1 α j ≥ 1 . A reasonable criterion for determining the optimal vaccine-allocation 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 k j=1 α j /E j , such that k j=1 α j ≥ 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 j-th vaccine per individual, for j = 1, . . . , k . Then a lower bound for the total cost of vaccination under the vaccination policy S is N(1 − 1 R 0 ) k j=1 α j c j /E j . Assuming N and 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 R 0 = 2.5 and six vaccines V 1 , . . . , 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 S V j represents a single-vaccine strategy such that only the j-th vaccine is used for the entire population; S n A , S n B , S n C , and S n D are multi-vaccine strategies based on a ranking-based method, such that n can take values 5, 10 or 20. For the rank-based 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 S , the total proportion of individuals who must be vaccinated under each strategy, relative effectiveness of strategies, ε(S) , and mean effectiveness of strategies, Ē S , respectively. The table body elements are α j for j = 1, . . . , 6. Table 1. Illustration of some strategies with their specific weights, herd-immunity threshold, relative effectiveness and mean effectiveness of strategies.  www.nature.com/scientificreports/ Table 1 is divided into two parts: the first represents single-variance strategies, while the second describes multi-vaccine strategies. The first part consists of six single-vaccine strategies such that p c j = P c S V j for j = 1, . . . , 6 . Therefore, p c 1 = 0.75 , p c 2 = 0.706 , p c 3 = 0.667 , p c 4 = 0.638 , p c 5 = 0.632 and p c 6 = 0.632 are critical vaccine thresholds for V 1 , . . . , V 6 , respectively. The least favourable strategy 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 multi-vaccine strategies can be proposed. However, only a few strategies are discussed, based on estimation methods described in the ranking-based method. There are four categories of multi-vaccine strategies in the second part of Table 1: • The first category of the second part is 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 α j = 1/6 ≈ 0.167 and p j = α j p c j , then p 1 ≈ 0.125 , p 2 ≈ 0.118 , p 3 ≈ 0.111 , p 4 ≈ 0.107 , p 5 ≈ 0.106 and p 6 ≈ 0.106 are proportions of the population that will be vaccinated by V 1 , . . . , V 6 to reach herd-immunity.
• 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 S 20 B proposes to use all vaccines, while S 20 C and S 20 D remove some vaccines from the vaccination program. • Category C of multi-vaccine 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 S 10 C and S 20 B 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, S 20 D assigns coefficients only to vaccines with largest effectiveness. • For a constant value of the shrinkage parameter n = n 0 , S n 0 D is better than S n 0 C and S n 0 C is better than S n 0 B , according to their mean and relative effectivenesses. When n is extremely large, all strategies in Categories B, C and D approach S 20 D , which is the most favourable strategy. • Note that relative effectiveness ε(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 multi-vaccine strategy to contain the infection.
Making decisions for vaccine allocation is a complex problem and depends on multi-factor priorities. Table 1 proposes plausible strategies that satisfy logistic, financial and social conditions of the public-health system. Overall, all strategies in Categories B, C, and D are better than S A , which is equivalent to a single-vaccine 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 π 1 , . . . , π m are the proportions of each variant, such that m j=1 π j = 1 . In addition, E 1 , . . . , E m are effectivenesses of the vaccine against variants. Under this setting, the vaccination reproduction number is R V =R 0 (1 − pĒ) , where R 0 = m j=1 π j R j 0 is the average reproduction number for all variants and Ē = m j=1 π j E j is the average effectiveness of the vaccine against the virus in the population. Note that estimating Ē 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 time-dependent; consequently, Ē 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 i-th vaccine against the j-th variant. Therefore, the vaccination reproduction number would be  www.nature.com/scientificreports/ Vaccination against symptomatic infection of COVID-19 in 2020-2021 has been deployed using BNT162b2 (Pfizer-BioNTech), mRNA-1273 (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 COVID-19 vaccine by July 31, 2021. This information is available on website of the government of Canada at https:// health-infob ase. 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. publi cheal thont ario. 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 Pfizer-BioNTech and Moderna, and between 60% and 70% for AstraZeneca. Since the effectiveness of vaccines decreases against VOC, we estimate the effectiveness of Pfizer-BioNTech, 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. publi shing. servi ce.
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., Ē i = j π j E 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:// outbr eak. info, which provides results of screening tests for mutations and whole-genome sequencing to assign COVID-19 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 single-dose and two-dose vaccination show "Vaccines not reported", which represent a group of people who are vaccinated by Pfizer-BioNTech, 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 single-dose and two-dose 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. www.nature.com/scientificreports/ • 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 Pfizer-BioNTech 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 R 0 = 2.5 for the original strain, we have R V = R 0 (1 − Partial p i E i1 − 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 R V =R 0 (1 − Partial p iĒi − Full p iĒi ) = 1.17 is our modified reproduction number. It follows that R V /R V ≈ 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 R V /R V may change substantially if the prevalence of new emerging vaccine-resistant 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 herd-immunity 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 vaccine-allocation strategy. First, we calculated the herd-immunity 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 single-vaccine strategy whose effectiveness is the sample mean of all effectivenesses, k j=1 E j /k , is not ideal, because many multi-vaccine strategies with smaller herd-immunity thresholds can be proposed. Third, our proposed modification of the herd-immunity threshold enables us to search for optimal combination of vaccines in order to reach herd immunity. We described methods to search for optimal vaccine-allocation strategies based on different priorities. Fourth, we showed that vaccine-allocation 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 herd-immunity 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 vaccine-acceptance rates for available vaccines for COVID-19 29,30 . Therefore, vaccine-acceptance rates or vaccine-hesitancy 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 time-dependent and can be modified as R V (t) = R 0 (1 − q(t) − k j=1 p j E 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 vaccine-induced 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 herd-immunity level. Therefore, a generalised condition for reaching herd immunity is k j=1 α j (Ē i /E i1 ) ≥ 1 − q(t)/(1 − 1/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 age-structured populations, for instance. The third is that they are not time-dependent. Further work is required to generalise our results in these circumstances. www.nature.com/scientificreports/ Vaccination at a level larger than the herd-immunity 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. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.