Vaccination and herd immunity thresholds in heterogeneous populations ===================================================================== * Elamin H. Elbasha * Abba B. Gumel ## Abstract It has been suggested, without rigorous mathematical analysis, that the classical vaccine-induced herd immunity threshold (HIT) assuming a homogeneous population can be substantially higher than the minimum HIT obtained when considering population heterogeneities. We investigated this claim by developing, and rigorously analyzing, a vaccination model that incorporates various forms of heterogeneity and compared it with a model of a homogeneous population. By employing a two-group vaccination model in heterogeneous populations, we theoretically established conditions under which heterogeneity leads to different HIT values, depending on the relative values of the contact rates for each group, the type of mixing between groups, relative vaccine efficacy, and the relative population size of each group. For example, under biased random mixing and when vaccinating a given group results in disproportionate prevention of higher transmission per capita, it is optimal to vaccinate that group before vaccinating other groups. We also found situations, under biased assortative mixing assumption, where it is optimal to vaccinate more than one group. We show that regardless of the form of mixing between groups, the HIT values assuming a heterogeneous population are always lower than the HIT values obtained from a corresponding model with a homogeneous population. Using realistic numerical examples and parametrization (e.g., assuming assortative mixing together with vaccine efficacy of 95% and basic reproduction number of 2.5), we demonstrate that the HIT value considering heterogeneity (e.g., biased assortative mixing) is significantly lower (40%) compared with a HIT value of (63%) assuming a homogeneous population. Keywords * Basic reproduction number * herd immunity threshold * heterogeneity * homogeneous population * mixing pattern * SVIR model ## 1. Introduction Since the identification of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and subsequent devastating impact of the 2019 novel coronavirus pandemic (COVID-19), a pertinent question being asked by the global public health and scientific research community has been what is the minimum fraction of the susceptible population that needs to be immunized for the pandemic to end? It has been widely reported that the herd immunity threshold (HIT) vaccination coverage, denoted by *v**, is given by the formula ![Graphic][1], where ℛ is *the basic reproduction number* of the model, defined as the average number of secondary infections caused by an infective person over his or her infectious period when introduced into a completely susceptible population [1]. For example, ℛ = 2.5 indicates the need to immunize more than 60% of the population to achieve HIT, so that the disease becomes eliminated from the population. However, this derivation for HIT is based on making several assumptions regarding the properties of the vaccine and characteristics of the population. Specifically, it assumes that the vaccine efficacy to protect against the acquisition of infection is perfect and lasts throughout the lifetime of the vaccinated individual (i.e., the vaccine does not wane). McLean and Blower [2], and other researchers [3,4], derived modified HIT formulae under various assumptions for vaccine properties. Their derivation is based on considering the following scenario for a cohort (childhood) vaccine. Suppose a hypothetical vaccine is efficacious in a fraction, *ε*, of recipients and confers full protection that wanes at a rate *ω per* unit time in a population with average life expectation of ![Graphic][2], so that the average duration of vaccine-induced protection is ![Graphic][3]. In this case, the critical vaccination coverage (i.e., the proportion of newborns that need to be vaccinated to achieve herd immunity, consequently leading to disease elimination), denoted by *v***, is given by [3] ![Formula][4] where *v** is as defined above. Thus, the herd immunity threshold (*v**) needs to be adjusted upward to reflect the imperfect efficacy of the vaccine (0 < *ε* < 1) and the fraction of the average lifetime a vaccinated individual remains protected (1/(*μ* + *ω*) ÷ 1/*μ*). For example, with ℛ= 2.5 and a vaccine with protective efficacy of only 70% (i.e., *ε* = 0.7) that lasts 90% of a lifetime (i.e., ![Graphic][5]) is used in the population, at least 95.2% of the population of the unvaccinated susceptible newborns needs to be vaccinated to achieve disease elimination (i.e., achieve HIT), since *v*** = 0.952 in this case. The above derivations assume a homogeneous, well-mixed population. However, the transmission of many infections (such as SARS-CoV-2) occurs in a diverse heterogeneous population. Hence, a more realistic approach to carry out the above computations will be to account for the relevant heterogeneities [5]. In other words, the computations need to be carried out for the case where the total population is divided into multiple groups with similar characteristics, such as age, contact patterns, infectious period, or social, cultural, demographic, or geographic factors. For example, several mathematical models for disease transmission employ different mixing patterns, such as those between different age groups [6,7]. These contact patterns are typically parametrized using empirical or synthetic social contact matrices estimated from population-based surveys [8]. ## 2. Formulation of Disease Transmission Model in Heterogeneous Populations The multigroup vaccination model for the transmission dynamics of a disease (such as SARS-CoV-2) in a heterogeneous population typically follows a standard susceptible-vaccinated-exposed-infected-recovered (SVEIR) Kermack-McKendrick-type compartmental modeling formulation [7, 9]. In this formulation, the total population at time t (denoted by *N*(*t*)) is sub-divided into *m* distinct homogeneous groups. Each group is further sub-divided into five disjoint (mutually exclusive) classes or compartments of unvaccinated susceptible (*S*(t)), vaccinated susceptible (*V*(*t*)), exposed (*E*(*t*)), infectious (*I*(*t*)), and recovered/removed (*R*(*t*)), so that *N**i* = *S**i* + *V**i* + *E**i* + *I**i* + *R**i*, with *i* = 1,2,…, *m*. Mathematically speaking, ‘exposed individuals’ are those who are newly infected with the pathogen but are not yet able to transmit the pathogen to other individuals (i.e., they are not infectious yet) [7]. The resulting SVEIR model, for the transmission dynamics of a disease in *m* heterogeneous groups or populations, is given by the following deterministic system of nonlinear differential equations (where a dot represents differentiation with respect to time *t*): ![Formula][6] ![Formula][7] where the *force of infection λ**i*(*t*) is given by ![Formula][8] with *β**i* as the transmission probability *per* contact for group *i, a**i* is the average number of contacts that an individual of group *i* has during a certain period of time (called group-specific activity level), and *c**ij* is the proportions of contacts that members of group *i* have with other individuals of group *j*. Mixing should meet the following closure relation [10]: ![Formula][9] That is, the total number of contacts that individuals of group *i* have with other individuals of group *j* during a certain period of time should be equal to the total number of contacts that individuals of group *j* have with other individuals of group *i*. In the model (1), heterogeneity between groups is captured through differences in demographic rates (i.e., birth and death rates), transmission probability per contact, contact rates, progression and recovery rates, and vaccine efficacy and waning rates. Adding all the equations of the model (1) gives the following equation for the rate of change of the total population: ![Formula][10] In the model (1), Λ*i* is the *per capita* recruitment (birth) into the population, *ω**i* is the vaccine waning rate, *λ**i*(*t*) is the force of infection, *μ**i* is the natural death rate (i.e., 1/*μ**i* is the average lifespan of a person in group *i*) and *ξ**i* is the vaccination rate. Furthermore, 0 < *ε**i* < 1 is the protective efficacy of the vaccine, *σ**i*is the rate at which exposed individuals develop clinical symptoms of the disease (i.e., 1/*σ**i* is the latency period of the disease), and γ*i* is the recovery rate. Some of the main assumptions made in the formulation of the model (1) are: 1. Within-group homogeneous mixing (i.e., although the model considers *m* heterogeneously mixed groups, contact patterns within each group is homogeneous). 2. Exponentially distributed waiting time in each epidemiological compartment. 3. The vaccine is not perfect (i.e., 0 < *ε**i* < 1), and the protection offered by the vaccine wanes over time (i.e., *ω* > 0). In addition, the vaccine has no therapeutic benefits. 4. No disease-induced mortality (so that the total population in each group remains constant). 5. Recovery induces permanent immunity against acquisition of future infection. ### 2.1 Disease-free equilibrium of general model in heterogeneous populations The model (1), for a single group, has been subjected to rigorous mathematical analysis in the literature. Specifically, results for its well-posedness, invariance of its solutions, and existence and asymptotic stability of equilibria (disease-free and endemic) have been established [7, 11]. The multigroup model (1) has a unique disease-free equilibrium given by ![Formula][11] It is convenient to assume that the population in each group *i* has reached a stationary (equilibrium) state such that ![Formula][12] It is also convenient to work with the fraction of the population in each group. The proportion of individuals in group *i* that are vaccinated (at the disease-free equilibrium) is given by: ![Formula][13] where, ![Formula][14] is the fraction of total individuals in group *i* relative to total population. ## 3. Analysis of Model in Heterogeneous Populations with Two Risk Groups Here, we focus on deriving analytic expressions for HITs for a special case of the heterogeneous population model (1) that considers only two risk groups (i.e., the model (1) with *m* = 2). For a two-group model, we first derive the vaccination reproduction number, ℛ*v*, and then find the minimum proportion of individuals in the community that need to be vaccinated in order to reduce ℛ*v* to a *v*alue less than 1, so that herd immunity is achieved. ### 3.1 Vaccination and basic reproduction numbers for two-group model The next-generation operator method [13] can be used to compute the vaccination reproduction number (and, subsequently, the basic reproduction number) of the special case of the model (1) with m = 2. It follows, using the notation in [12], that the non-negative matrix of new infection terms (*F*) and the M-Matrix (*V*) of linear transition terms in the infected compartments are given, respectively, by (where ![Graphic][15] and ![Graphic][16] are the total population of group 1 and group 2, respectively; similarly, ![Graphic][17] and ![Graphic][18] represent the HIT for groups 1 and 2, respectively): ![Formula][19] from which it follows that the vaccination reproduction number is given by (where *ρ* is the spectral radius) ![Formula][20] where, ![Formula][21] ![Formula][22] In deriving equation (3), we utilized the following definition of the constituent basic reproduction numbers associated with disease transmission between individuals in group *i* with individuals in group *j*: ![Formula][23] where the index *i* represents group *i* and the index *j* represents group *j*. To obtain the basic reproduction number (ℛ) associated with the two-group model, we set the vaccination coverage rates in the expression for ℛ*v*, given by (3), to zero (i.e., we set ![Graphic][24] in (3)). This gives, ![Formula][25] ### 3.2 Herd immunity thresholds for a two-group model with heterogeneous populations For herd immunity threshold determination or computation of a two-group model, the objective is to find the values of the respective HITs, ![Graphic][26] and ![Graphic][27], such that the total vaccine coverage (i.e., the proportion of individuals in the community that is vaccinated), given by ![Formula][28] is at its minimum and vaccinated reproduction number, ℛ*v*, given by (3), is less than or equal to one. Formally, the optimization problem can be written as choosing ![Graphic][29] and ![Graphic][30] to ![Formula][31] subject to, ![Formula][32] where ℛ*v* is given by equation (3). The solution of this nonlinear optimization problem will be characterized using a geometrical approach. Specifically, we compare the shape of the curve depicting values of vaccination coverage where the vaccinated reproduction number is equal to one (ℛ*v* = 1) with the contour lines (or level sets) ![Graphic][33] (as illustrated in Figure 1). Each contour line represents the locus of vaccination coverage combinations ![Graphic][34] that yield the same level of total vaccination coverage at the population level. The blue contour lines correspond to lower and lower levels of total vaccination coverage when moving in the southwestern direction toward the origin. ![](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/03/01/2021.02.26.21252553/F1/graphic-20.medium.gif) [](http://medrxiv.org/content/early/2021/03/01/2021.02.26.21252553/F1/graphic-20) ![](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/03/01/2021.02.26.21252553/F1/graphic-21.medium.gif) [](http://medrxiv.org/content/early/2021/03/01/2021.02.26.21252553/F1/graphic-21) Figure 1. Vaccination threshold values for group 1 and group 2. The orange curve shows values of vaccination coverage where the vaccinated reproductive number is equal to one: ℛ*v* = 1. The slope of this curve when intersects with the *y-*axis is given by ![Graphic][35]. The blue level curves show different values of total vaccination coverage going down in the direction of the origin: ![Graphic][36]. The slope of these level curves is – ![Graphic][37]. Under biased random mixing (a) when the blue line is flatter than the orange curve when it intersects the *y-*axis, the closest blue line to the origin that intersect the orange curve happens when ![Graphic][38]: (b) when the blue line is steeper than the orange curve when it intersects the *x-* axis, the closest blue line to the origin that intersect the orange curve happens when ![Graphic][39]. Under biased assortative mixing (c) the optimum occurs when the blue line is tangent to the orange curve. In the homogeneous population model ![Graphic][40], the optimum occurs at the intersection of the black 45-degree line with the orange curve. The solution of the equation of the orange curve ℛ*v* = 1 yields a value of ![Graphic][41] as a function of ![Graphic][42]: ![Formula][43] It follows from equation (6) that this function intersects the *x-*axis ![Graphic][44] and *y-*axis ![Graphic][45] at ![Formula][46] respectively. The slope of the orange curve (Figure 1) is given by ![Formula][47] The slope of the orange curve, at the point of intersection with the *y-*axis (i.e., ![Graphic][48]), is given by ![Formula][49] and the slope of the orange curve when it intersects with the *x-*axis (i.e., ![Graphic][50]) is given by ![Formula][51] The equations of the blue level curves are: ![Graphic][52], where *A* is an arbitrary constant (Figure 1). The slope of these level curves is ![Formula][53] Depending on the situation, the geometrical approach involves comparing the slopes of the curves given by equations (8)–(11). We consider the following two scenarios, based on which group(s) the vaccine is prioritized to. ### Scenario 1: Vaccinating only one group We start with the case where ℛ*v* can be brought down to one by vaccinating either group 1 alone or group 2 alone. That is, the values of ![Graphic][54] and ![Graphic][55] given in equation (7) are less than 1. We note that the type of mixing between groups (i.e., *c**ij*) determines the sign of the expression ℛ12ℛ21 − ℛ11ℛ22. By using equation (4), we can establish that ![Formula][56] We distinguish two types of mixing. #### A) Biased random mixing: ℛ12ℛ 21 ≥ ℛ11ℛ 22 When transmission occurs more because of mixing between groups rather than because of mixing within the groups, we call this type of mixing *biased* toward *random mixing*. The *separable proportionate mixing* is a special type of biased random mixing [7, 10]. i) Vaccinating group 1 disproportionately contributes more to prevention of *per-capita* transmission. Here, the following inequality holds: ![Formula][57] Figure 1a depicts the situation where the absolute value of the slope of the orange curve when it intersects the *y*-axis (equation (9)) is greater than the absolute value of the slope of the blue line (equation (11)) (i.e., blue line is flatter). That is (hence the inequality above), ![Formula][58] In this case, the closest blue line to the origin that intersect the orange curve happens when ![Formula][59] and the herd immunity threshold is given by ![Formula][60] (ii) Vaccinating group 1 disproportionately contributes less to prevention of per-capita transmission: ![Formula][61] In this case, the absolute value of the slope of the orange curve when it intersects the y-axis (equation (9)) is less than or equal the absolute value of the slope of the blue line (equation (11)) (i.e., blue line is steeper) (illustrated in Figure 1b), ![Formula][62] and the herd immunity threshold is ![Formula][63] In summary, the above results show that, under the biased random mixing assumption, achieving herd immunity entails exclusively optimizing vaccination coverage among the group that results in relatively more prevention of *per-capita* disease transmission. #### B) Biased assortative mixing: ℛ12ℛ 21 < ℛ11ℛ 22 When transmission occurs more due to mixing within groups rather than due to mixing between groups, we refer to this type of mixing as *biased* toward *assortative mixing*. In this case, the optimization problem has both boundary (corner) and interior solutions. We start with the two boundary solutions followed by the interior solution. i) Vaccinating group 1 disproportionately contributes more to prevention of *per-capita* transmission: ![Formula][64] In this case, the absolute values of the slope of the orange curve when it intersects the *y-*axis (equation (9)) and *x-*axis (equation (10)) are greater than the absolute value of the slope of the blue line (equation (11)) (i.e., blue line is flatter). In this case, ![Formula][65] and the herd immunity threshold is given by equation (12). In other words, we achieve optimal results by allocating all necessary vaccine resources to group 1. ii) Vaccinating group 1 disproportionately contributes less to prevention of *per-capita* transmission. Here, we have: ![Formula][66] In this case the absolute values of the slope of the orange curve when it intersects the *y-*axis (equation (9)) and *x-*axis (equation (10)) are less than the absolute value of the slope of the blue line (equation (11)) (i.e., blue line is steeper). In this case, ![Formula][67] and the herd immunity threshold is given by equation (13). Here, optimal results are achieved by allocating all vaccine resources to group 2. ### Scenario 2: Vaccinating both groups This scenario includes one interior solution and 2 boundary solutions. For this scenario, we do not need to assume that ℛ*v* can be brought down to one by vaccinating either group 1 alone or group 2 alone. i) Interior solution. For this solution to occur, the following inequalities must hold: ![Formula][68] When the absolute value of the slope of the blue line (equation (11)) is between the absolute values of the slope of the orange curve when it intersects the *y-*axis (equation (9)) and *x-*axis (equation (10)), it is optimal to vaccinate both groups (illustrated in Figure 1c). The interior solution is obtained when the slope of the orange curve (equation (8)) is equal to the slope of the blue curve (equation (11)). Thus, by equating the right-hand side of equations (8) and (11) and solving for ![Graphic][69] and using the resulting solution in equations (6), we have: ![Formula][70] It follows from the assumption of biased assortative mixing and the inequalities above a positive **interior** solution exists if ![Formula][71] The herd immunity threshold is ![Formula][72] ii) Boundary solutions. The case where ℛ*v* cannot be brought down to one by vaccinating either group 1 alone or group 2 alone (i.e., the values of either ![Graphic][73] or ![Graphic][74] are greater than 1) always leads to boundary solutions. For example, when ![Graphic][75] and ![Graphic][76] and ![Graphic][77], we first maximize vaccination coverage among group 2 and then find the value of coverage among group 1 that brings ℛ*v* to one. In this case, ![Formula][78] and the herd immunity threshold is ![Formula][79] Similarly, when![Graphic][80], we may have, ![Formula][81] and the herd immunity threshold is ![Formula][82] Thus, the optimal solution and HIT are summarized as follows: 1. an interior solution if inequalities given in (14) are satisfied. The corresponding HIT value is given by equation (15); 2. four boundary solutions exist if any of inequalities given in (14) does not hold. The associated HIT value is given by equations (12), (13), (16), and (17); 3. no solution otherwise. These analyses show that, for a two-group vaccination model in heterogeneous populations, the optimum vaccination program depends on the relative values of the constituent reproduction numbers, ℛ12, ℛ21,ℛ11, ℛ22, relative vaccine efficacy, and the relative population sizes of the two groups. The values of the constituent reproduction numbers are determined by the type of mixing allowed or assumed between the two groups. When the mixing between the groups is biased towards random mixing, achieving herd immunity entails restricting vaccination coverage to the group that results in relatively more prevention of *per-capita* transmission. If herd immunity cannot be achieved by vaccinating one of the two groups alone, vaccination coverage among the group that results in relatively more prevention of *per-capita* transmission should be maximized before vaccinating the other group. These scenarios occur under biased assortative mixing too, but scenarios involving vaccinating both groups are more common, including interior optimum where coverage of any of the two groups is less than 100%. ## 4. Analysis of Model with a Homogeneous Population The model (1) can be reduced to that with homogeneous population by assuming that each individual in any of the two groups is identical with every other individual in the community. That is, we achieve homogeneity by setting *c**ij*= *c, a**i* = *a, ξ**i* = *ξ, ε**i* = *ε, μ**i* = *μ*, ⋀*i* = ⋀, *σ**i* = *σ*, γ*i* = γ for all *i* and *j* into the model (1). Specifically, the vaccination model with a homogeneous population is obtained from system (1) by dropping the group subscript *i* and re-defining the force of infection as ![Formula][83] Using the next-generation operator method [13], the *vaccination reproduction number* associated with the resulting homogeneous model, denoted by ℛ*v*, is given by: ![Formula][84] where the basic reproduction number, ℛ, is given by ![Formula][85] Since the vaccination coverage at the disease-free equilibrium (*v**) for the multigroup model with homogeneous population is ![Formula][86] it follows, by using (18) in (20), that ![Formula][87] It should be noted that this relationship between ℛ*v* and ℛ (equation (21)) in the homogeneous population model can obtained directly from the formula of the vaccination reproduction number of the heterogeneous population model by substituting *ε**i* = *ε* and ![Graphic][88], *i* = 1,2, into equation (3) and using the definition of ℛ given by equation (5). Setting ℛ*v* to one and solving for *v** gives the following threshold value needed to achieve herd immunity for the model with homogeneous population [3, 7]: ![Formula][89] Thus, in constructing the homogeneous population model corresponding to the heterogeneous model (1) for comparison, we followed previous studies and matched the two modeling types (i.e., homogeneous vs. heterogeneous population) according to the expressions or values of their respective basic reproduction number, ℛ [14,15]. ### 4.1 Comparison of HIT values using homogeneous and heterogeneous population models under proportionate mixing One of the simplest types of mixing in disease epidemiology is the separable proportionate mixing, in which the contacts of a person of group *i* are distributed over those of other groups in proportion to the activity levels and sizes of the other groups [7]. Thus, with proportionate mixing, proportions of contacts that members of group *i* have with group *j, c**ij*, is given by ![Formula][90] where *a**i* and ![Graphic][91] are as defined before. By substituting this definition of *c**ij* into the formulae for the constituent reproduction numbers in equation (4), it can be shown that the assumption of proportionate mixing implies that ℛ12ℛ21 = ℛ11ℛ22, which in turns implies, Δ2= 0. Thus, for this scenario of proportionate mixing, the vaccination reproduction number reduces to, ![Formula][92] We now seek to answer the question: what are the values of ![Graphic][93] and ![Graphic][94] such that total vaccination coverage, ![Graphic][95], is at its minimum and *R**v* ≤ 1? The answer depends on the relationship between the ratio of constituent reproduction numbers ![Graphic][96] adjusted by efficacy and the ratio ofpopulation ![Graphic][97] in the two groups. Solution of this simple linear programming problem is a special case of the biased random mixing where ℛ12ℛ21 = ℛ11ℛ22. As before, there are three scenarios: **Scenario (i)** Vaccinating group 1 disproportionately contributes more to prevention of *per-capita* transmission: ![Graphic][98]. In this case, ![Formula][99] and the herd immunity threshold is ![Formula][100] **Scenario (ii)** Vaccinating group 1 disproportionately contributes less to prevention of *per-capita* transmission: ![Graphic][101]. In this case, ![Formula][102] and the herd immunity threshold is ![Formula][103] **Scenario (iii)** Vaccinating both groups contribute equally to prevention of transmission: ![Graphic][104] In this case, values of ![Graphic][105] and ![Graphic][106] such that ![Formula][107] will yield the minimum fraction that need to be vaccinated to achieve herd immunity. To facilitate the comparison of the herd immunity thresholds between homogeneous and heterogeneous population models, we need to make sure efficacy in the two models is the same. One approach is to assume that efficacy does not vary across the two groups such that *ε*1 = *ε*2 = *ε*. If ![Graphic][108], the threshold vaccine coverage under heterogeneous population model is (given by the right-hand side of equation (23)): ![Formula][109] and that for the homogeneous population model is given by ![Formula][110] since ℛ = ℛ11 + ℛ22 > 1 under proportionate mixing. It can be shown that ![Formula][111] Upon simplifications, the above inequality holds if ![Formula][112] Noting (and using) our starting assumption ![Graphic][113], it follows that ![Graphic][114] or ![Formula][115] Therefore, it follows from the above inequality, that the threshold vaccine coverage in the heterogeneous population model, under scenario (i), is always less than that in the corresponding homogeneous population model. Similarly, if ![Graphic][116], we can follow the same approach above and show that the threshold vaccine coverage under heterogeneous population model, under scenario (ii), given by equation (24) is always lower than that under the homogeneous population: ![Formula][117] When ![Graphic][118], it follows, under equal efficacy, that ![Graphic][119] or ![Formula][120] Thus, the threshold vaccine coverage under heterogeneous population model, under scenario (iii), given by equation (25) is always equal to that under the homogeneous population: ![Formula][121] Therefore, under the form of proportionate mixing between groups, we show analytically that the HIT value in the heterogeneous population model is always less than or equal to the HIT value in a homogeneous population model. ### 4.2 Comparison of HIT values using homogeneous and heterogeneous population models under general mixing To show that the HIT value in the heterogeneous population model is always less than or equal to the HIT value in a homogeneous population model under general mixing between groups we utilize the geometrical approach depicted in Figure 1. Recall that the blue level curve represents total vaccination coverage. Thus, line *AA* corresponds to HIT value in heterogeneous population model whereas *BB* corresponds to HIT value in the homogeneous population model. We consider scenarios with two boundary solutions and one interior solution: 1. Vaccinating group 1 only (Figure 1a). In this case, line *AA* is closer to the origin than line *BB*. Hence, the HIT value in heterogeneous population model is lower than the HIT value in the homogeneous population model. 2. Vaccinating group 2 only (Figure 1b). In this case, line *AA* is closer to the origin than line *BB*. Hence, the HIT value in heterogeneous population model is lower than the HIT value in the homogeneous population model. 3. Vaccinating both group 1 and group 2 (Figure 1c). In this case, line *AA* is closer to the origin than line *BB*. Hence, the HIT value in heterogeneous population model is lower than the HIT value in the homogeneous population model. ## 5. Numerical Analysis of HIT Values in Heterogeneous and Homogeneous Population Models Figure 2 illustrates numerically the different scenarios leading to different HIT values in the heterogeneous population model and compare them with the HIT values in a corresponding homogeneous population model. The orange curve shows values of vaccination coverage where ℛ*v* = 1 and the blue level curves show different values of total vaccination coverage going down in the direction of the origin ![Graphic][122]. In a homogeneous population model, vaccination coverage (determined by the intersection with the orange curve) is uniform across the two groups as shown by the dotted black 45-degree line. In all scenarios we set vaccine efficacy to 95% across the two groups. ![Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/03/01/2021.02.26.21252553/F2.medium.gif) [Figure 2.](http://medrxiv.org/content/early/2021/03/01/2021.02.26.21252553/F2) Figure 2. Illustration of optimal least vaccine coverage determination by group that satisfies the constraint ℛ*v* = 1 (orange curve). Parameter values: *ε* = 0.95, *ε*1 = 0.95, *ε*2 = 0.95. (a–c): ℛ12 = 1.0, ℛ21 = 0.8, ℛ22 = 1.3, ![Graphic][123]. a) ℛ11 = 0.5, ℛ = 1.88; b) ℛ11 = 2, ℛ = 2.6; c) ℛ11 = 1.2, ℛ = 2.16; d) ℛ12 = 0.8, ℛ21 = 1.0, ℛ22 = 0.5, ![Graphic][124]; e)ℛ12 = 0.8, ℛ21 = 1.0, ℛ22 = 0.5, ![Graphic][125]. 1. The chosen parameters values represent a situation of biased random mixing between groups, and ![Graphic][126] and ![Graphic][127]. As a result, the blue line is steeper than the orange curve when the latter intersects the *y-*axis, and the closest blue line to the origin that intersect the orange curve happens when ![Graphic][128] (Figure 2a). Given that group 2 represents only 25% of the population, the overall HIT value is just 17.2% compared with the HIT value in a homogeneous population model of 49.3%. 2. The chosen parameters values represent a situation of biased assortative mixing between groups, ![Graphic][129] and ![Graphic][130], and conditions (14) is satisfied. As a result, it is optimal to vaccinate both groups ![Graphic][131] for an overall HIT value of 64.3% (Figure 2b). The HIT value in a homogeneous population model is 64.9%. 3. The chosen parameters values represent a situation of biased assortative mixing between groups, and ![Graphic][132] and ![Graphic][133]. As a result, the blue line is steeper than the orange curve, and it is optimal to vaccinate all of group 2 and 20.6% of group ![Graphic][134] for an overall HIT value of 40.4% (Figure 2c). The HIT value in a homogeneous population model is 56.2%. 4. The chosen parameters values represent a situation of biased random mixing between groups and ![Graphic][135] and ![Graphic][136]. As a result, the blue line is flatter than the orange curve, and it is optimal to vaccinate group 1 only ![Graphic][137] for an overall HIT value of 13.5% (Figure 2d). The HIT value in a homogeneous population model is 47.1%. 5. The chosen parameters values represent a situation of biased assortative mixing between groups and ![Graphic][138] and ![Graphic][139]. As a result, the blue line is flatter than the orange curve, and it is optimal to vaccinate all of group 1 and 20.6% of group ![Graphic][140] for an overall HIT value of 40.5% (Figure 2e). The HIT value in a homogeneous population model is 63.3%. Figure 3 illustrates different optimal solutions and HIT values for various values of the basic reproduction numbers. The figure shows the situation where relative to its small size group 2 is contributing more to transmission for low values of ℛ11and ℛ and there is a need to vaccinate more people from group 2. As ℛ11(and ℛ) increases, necessary vaccination coverage among group 2 increases until all of group 2 is vaccinated. As ℛ11(and ℛ) increases further, both groups are vaccinated, but vaccination coverage among group 1 increases whereas that among group 2 decreases. Of note, the herd immunity threshold for the homogeneous population (red curve) is consistently higher, but the difference between the two shrinks as the basic reproduction number increases. ![Figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2021/03/01/2021.02.26.21252553/F3.medium.gif) [Figure 3.](http://medrxiv.org/content/early/2021/03/01/2021.02.26.21252553/F3) Figure 3. Vaccine coverage above which herd immunity is achieved by group and basic reproduction number. Parameters values: ℛ12 = 1.0, ℛ21 = 0.8, ℛ22 = 1.3, *ε* = 0.95, *ε*1 = 0.95, *ε*2 = 0.95, ![Graphic][141]. Using the definition of ℛ, we chose ![Graphic][142]. ## 6. Discussion Our theoretical results, for a two-group model, suggest that deriving the vaccine-induced herd immunity threshold in the heterogeneous population model is much more complicated than for the model with homogeneous population. The herd immunity threshold for each vaccinated group depends on the relative values of the constituent reproduction numbers, relative vaccine efficacy, and the relative population sizes of the two groups. The values of the reproduction numbers are determined by the level and duration of infectiousness of a contact for each group, contact rates for each group, as well as the type of mixing between groups. We show that, under biased random mixing and when vaccinating a given group results in disproportionately prevention of higher transmission *per capita*, it is optimal to optimize vaccination of that group before vaccinating other groups. We also found situations, under biased assortative mixing assumption, where it is optimal to vaccinate more than one group. We show that population heterogeneities tend to result in lower HIT values, compared with the homogeneous population case. This is true under both proportionate and other types of mixing among heterogeneous populations. Using realistic numerical examples and parametrization (e.g., assuming biased assortative mixing with vaccine efficacy of 95% and basic reproduction number of 2.5) we illustrate this finding, where the HIT value considering heterogeneity is shown to be significantly lower (40%) compared with a HIT value assuming a homogeneous population (63%). Although our analysis utilized a two-group model, our findings can be extended to models of multiple groups. Although more complicated than a two-group model, rigorous analyses of more realistic models with many heterogeneous groups can be conducted using the methods included this paper. ## Data Availability The work is theoretical. No new data has been generated. ## Acknowledgements ABG acknowledges the support, in part, of the Simons Foundation (Award #585022) and the National Science Foundation (Award #1917512). * Received February 26, 2021. * Revision received February 26, 2021. * Accepted March 1, 2021. * © 2021, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NoDerivs 4.0 International), CC BY-ND 4.0, as described at [http://creativecommons.org/licenses/by-nd/4.0/](http://creativecommons.org/licenses/by-nd/4.0/) ## References 1. 1.Randolph H, Barreiro L (2020) Herd Immunity: Understanding COVID-19. Immunity 52(5):737–741 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.immuni.2020.04.012&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32433946&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F03%2F01%2F2021.02.26.21252553.atom) 2. 2.McLean AR, Blower SM (1993) Imperfect vaccines and herd immunity to HIV. Proc R Soc Lond B 253:9–13. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1098/rspb.1993.0075&link_type=DOI) 3. 3.Scherer A, McLean A (2002) Mathematical models of vaccination. Br Med Bull 62:187–99 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bmb/62.1.187&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12176860&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F03%2F01%2F2021.02.26.21252553.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000177294400015&link_type=ISI) 4. 4.Magpantay FMG (2017) Vaccine impact in homogeneous and age-structured models. J Math Biol 75(6-7): 1591–1617 5. 5.Britton T, Ball F, Trapman P (2020) A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2. Science 369(6505):846–849 [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzNjkvNjUwNS84NDYiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMS8wMy8wMS8yMDIxLjAyLjI2LjIxMjUyNTUzLmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 6. 6.1. D. Mollison, ed Jacquez JA, Simon CP, and Koopman JS (1996) Core groups and the R0’s for subgroups in heterogeneous SIS models, in Epidemic Models: Their Structure and Relation to Data, D. Mollison, ed., Cambridge University Press, Cambridge, UK, pp. 279–301 7. 7.Hethcote HW (2000). The mathematics of infectious diseases. SIAM Rev 42:599–653 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1137/S0036144500371907&link_type=DOI) 8. 8.Mossong J, Hens N, Jit M, Beutels P, Auranen K, Mikolajczyk R, et al. (2008) Social Contacts and Mixing Patterns Relevant to the Spread of Infectious Diseases. PLoS Med 5:e74 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pmed.0050074&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18366252&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F03%2F01%2F2021.02.26.21252553.atom) 9. 9.Kermack WO, McKendrick AG (1927) Contributions to the mathematical theory of epidemics, part 1, Proc. Roy. Soc. London Ser. A, 115, pp. 700–721 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1098/rspa.1927.0118&link_type=DOI) 10. 10.Glasser J, Feng Z, Moylan A, Del Valle S, Castillo-Chavez C (2012) Mixing in age-structured population models of infectious diseases. Math Biosci 235(1):1–7 [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22037144&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F03%2F01%2F2021.02.26.21252553.atom) 11. 11.Gumel AB, Iboi EA, Ngonghala CN, Elbasha EH (2020) A primer on using mathematics to understand COVID-19 dynamics: Modeling, analysis and simulations. Infect Dis Model 6:148–168 12. 12.van den Driessche P, Watmough J (2002) Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci 180:29–48 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0025-5564(02)00108-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12387915&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F03%2F01%2F2021.02.26.21252553.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000179220600004&link_type=ISI) 13. 13.Diekmann O, Heesterbeek JA, Roberts MG (2010) The construction of next-generation matrices for compartmental epidemic models. J R Soc Interface 7(47):873–885 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1098/rsif.2009.0386&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19892718&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F03%2F01%2F2021.02.26.21252553.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000276992100002&link_type=ISI) 14. 14.Andreasen V (2011) The final size of an epidemic and its relation to the basic reproduction number. Bull Math Biol 73:2305–2321 [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s11538-010-9623-3&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21210241&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2021%2F03%2F01%2F2021.02.26.21252553.atom) 15. 15.Clancy D, Pearce Christopher J (2013) The effect of population heterogeneities upon spread of infection. J Math Biol 67:963–987 [1]: /embed/inline-graphic-1.gif [2]: /embed/inline-graphic-2.gif [3]: /embed/inline-graphic-3.gif [4]: /embed/graphic-1.gif [5]: /embed/inline-graphic-4.gif [6]: /embed/graphic-2.gif [7]: /embed/graphic-3.gif [8]: /embed/graphic-4.gif [9]: /embed/graphic-5.gif [10]: /embed/graphic-6.gif [11]: /embed/graphic-7.gif [12]: /embed/graphic-8.gif [13]: /embed/graphic-9.gif [14]: /embed/graphic-10.gif [15]: /embed/inline-graphic-5.gif [16]: /embed/inline-graphic-6.gif [17]: /embed/inline-graphic-7.gif [18]: /embed/inline-graphic-8.gif [19]: /embed/graphic-11.gif [20]: /embed/graphic-12.gif [21]: /embed/graphic-13.gif [22]: /embed/graphic-14.gif [23]: /embed/graphic-15.gif [24]: /embed/inline-graphic-9.gif [25]: /embed/graphic-16.gif [26]: /embed/inline-graphic-10.gif [27]: /embed/inline-graphic-11.gif [28]: /embed/graphic-17.gif [29]: /embed/inline-graphic-12.gif [30]: /embed/inline-graphic-13.gif [31]: /embed/graphic-18.gif [32]: /embed/graphic-19.gif [33]: /embed/inline-graphic-14.gif [34]: /embed/inline-graphic-15.gif [35]: F1/embed/inline-graphic-16.gif [36]: F1/embed/inline-graphic-17.gif [37]: F1/embed/inline-graphic-18.gif [38]: F1/embed/inline-graphic-19.gif [39]: F1/embed/inline-graphic-20.gif [40]: F1/embed/inline-graphic-21.gif [41]: /embed/inline-graphic-22.gif [42]: /embed/inline-graphic-23.gif [43]: /embed/graphic-22.gif [44]: /embed/inline-graphic-24.gif [45]: /embed/inline-graphic-25.gif [46]: /embed/graphic-23.gif [47]: /embed/graphic-24.gif [48]: /embed/inline-graphic-26.gif [49]: /embed/graphic-25.gif [50]: /embed/inline-graphic-27.gif [51]: /embed/graphic-26.gif [52]: /embed/inline-graphic-28.gif [53]: /embed/graphic-27.gif [54]: /embed/inline-graphic-29.gif [55]: /embed/inline-graphic-30.gif [56]: /embed/graphic-28.gif [57]: /embed/graphic-29.gif [58]: /embed/graphic-30.gif [59]: /embed/graphic-31.gif [60]: /embed/graphic-32.gif [61]: /embed/graphic-33.gif [62]: /embed/graphic-34.gif [63]: /embed/graphic-35.gif [64]: /embed/graphic-36.gif [65]: /embed/graphic-37.gif [66]: /embed/graphic-38.gif [67]: /embed/graphic-39.gif [68]: /embed/graphic-40.gif [69]: /embed/inline-graphic-31.gif [70]: /embed/graphic-41.gif [71]: /embed/graphic-42.gif [72]: /embed/graphic-43.gif [73]: /embed/inline-graphic-32.gif [74]: /embed/inline-graphic-33.gif [75]: /embed/inline-graphic-34.gif [76]: /embed/inline-graphic-35.gif [77]: /embed/inline-graphic-36.gif [78]: /embed/graphic-44.gif [79]: /embed/graphic-45.gif [80]: /embed/inline-graphic-37.gif [81]: /embed/graphic-46.gif [82]: /embed/graphic-47.gif [83]: /embed/graphic-48.gif [84]: /embed/graphic-49.gif [85]: /embed/graphic-50.gif [86]: /embed/graphic-51.gif [87]: /embed/graphic-52.gif [88]: /embed/inline-graphic-38.gif [89]: /embed/graphic-53.gif [90]: /embed/graphic-54.gif [91]: /embed/inline-graphic-39.gif [92]: /embed/graphic-55.gif [93]: /embed/inline-graphic-40.gif [94]: /embed/inline-graphic-41.gif [95]: /embed/inline-graphic-42.gif [96]: /embed/inline-graphic-43.gif [97]: /embed/inline-graphic-44.gif [98]: /embed/inline-graphic-45.gif [99]: /embed/graphic-56.gif [100]: /embed/graphic-57.gif [101]: /embed/inline-graphic-46.gif [102]: /embed/graphic-58.gif [103]: /embed/graphic-59.gif [104]: /embed/inline-graphic-47.gif [105]: /embed/inline-graphic-48.gif [106]: /embed/inline-graphic-49.gif [107]: /embed/graphic-60.gif [108]: /embed/inline-graphic-50.gif [109]: /embed/graphic-61.gif [110]: /embed/graphic-62.gif [111]: /embed/graphic-63.gif [112]: /embed/graphic-64.gif [113]: /embed/inline-graphic-51.gif [114]: /embed/inline-graphic-52.gif [115]: /embed/graphic-65.gif [116]: /embed/inline-graphic-53.gif [117]: /embed/graphic-66.gif [118]: /embed/inline-graphic-54.gif [119]: /embed/inline-graphic-55.gif [120]: /embed/graphic-67.gif [121]: /embed/graphic-68.gif [122]: /embed/inline-graphic-56.gif [123]: F2/embed/inline-graphic-57.gif [124]: F2/embed/inline-graphic-58.gif [125]: F2/embed/inline-graphic-59.gif [126]: /embed/inline-graphic-60.gif [127]: /embed/inline-graphic-61.gif [128]: /embed/inline-graphic-62.gif [129]: /embed/inline-graphic-63.gif [130]: /embed/inline-graphic-64.gif [131]: /embed/inline-graphic-65.gif [132]: /embed/inline-graphic-66.gif [133]: /embed/inline-graphic-67.gif [134]: /embed/inline-graphic-68.gif [135]: /embed/inline-graphic-69.gif [136]: /embed/inline-graphic-70.gif [137]: /embed/inline-graphic-71.gif [138]: /embed/inline-graphic-72.gif [139]: /embed/inline-graphic-73.gif [140]: /embed/inline-graphic-74.gif [141]: F3/embed/inline-graphic-75.gif [142]: F3/embed/inline-graphic-76.gif