When mathematical models are used to evaluate vaccination strategies, results have been sensitive to population heterogeneity. [1][2][3][4] A quantity that is commonly used to design disease control policies is the basic reproduction number, denoted by R 0 , which has the property that outbreaks are possible if and only if R 0 > 1. If programs for disease control are implemented, the corresponding quantity is the effective reproduction number, denoted by R e . Typically R e < R 0 , as we consider only effective disease control strategies.
The meta-population model presented in this section was studied by Feng et al., 1 who demonstrated how population heterogeneity may affect vaccination strategies. This is an susceptible-infected-recovered (SIR) type of compartmental model for the transmission of pathogens causing infectious diseases when immunity by virtue of either infection or immunization is permanent (but the same framework can be applied to other models including SEIR, SIS, etc.). The meta-population (a population of populations) consists of n subgroups or sub-populations connected by a mixing matrix or network, which describes person-to-person contacts within and between subgroups.
Let S i (t), I i (t), R i (t) denote the numbers of susceptible, infectious, recovered/immune individuals in group i at time t. Then N i = S i + I i + R i is the total population size of group i, and p i represents the proportion of group i that is immunized, i = 1, 2, . . . , n. The model reads:
with initial conditions S i (0) > 0, I i (0) ≥ 0 and > 0 for some i, R i (0) ≥ 0. The parameter γ is the per capita recovery rate, µ is the per capita birth and death rates, which are assumed to be equal (so that the total population in group i remains constant). The function λ i (t) represents the force of infection for susceptible members of group i given by
where a i denotes the average number of contacts an individual in sub-population i has per unit of time (which represents the activity level of group i), and σ is the probability of infection per contact with an infectious individual. The parameter c ij is the proportion of the contacts that members of the ith sub-population have with members of the jth, defined by
where i is the fraction of their contacts that members of the ith sub-population reserve for others in their own sub-population (called preference), δ ij is the Kronecker data,
is the contribution of members of the jth sub-population to the pool of unreserved contacts, and a j is the per capita contact rate (called activity). The function c ij allows mixing to range from proportionate (all i = 0) to preferential (some 0 < i ≤ 1). 1,2,[4][5][6][7] A fundamental difference between models with and without heterogeneity in terms of optimal immunization strategy is illustrated in Feng et al. 1 For example, consider the case of n = 2 groups. In the absence of heterogeneity (i.e., parameters for the two groups are identical), p i = p (i = 1, 2), and the basic and effective reproduction numbers for both groups and the meta-population are
The threshold population immunity is
such that R e < 1 if and only if p > p c , which eventually leads to disease elimination. That is, for disease elimination, the optimal immunization proportions in the two sub-populations are p i > p ic with p ic being given in (1.6).
When heterogeneity is considered, however, p 1c and p 2c can differ. In fact, the basic reproduction number for group i is R 0i = σa i /(γ + µ), and the effective reproduction number for group i is
But, for the metapopulation, the effective reproduction number is given by
where
(1.8)
Any combination of vaccine coverage (p 1 , p 2 ) for which R e (p 1 , p 2 ) < 1 can lead to elimination of disease. Feng et al. 1 showed that there are infinitely many critical combinations (p 1c , p 2c ) that have the property that R e (p 1 , p 2 ) < 1 for all p 1 > p 1c and p 2 > p 2c , including the homogeneous allocation
The minimum number of vaccine doses required to achieve disease elimination depends on the mixing matrix (c ij ). Most interestingly, the homogeneous allocation given in (1.9) may not be the optimal strategy. The thick (red) curve in Fig. 1 corresponds to R e (p 1c , p 2c ) = 1. This is for the case of N 1 = N 2 , so the total vaccine doses p 1 N 1 + p 2 N 2 = (p 1 + p 2 )N 1 is minimized if and only if p 1 + p 2 is minimized. All points (p 1 , p 2 ) on the dotted line (tangent to the R e = 1 curve at the star) correspond to same doses, including the homogeneous allocation (p * 1c , p * 2c ) (intersection of the dotted and dashed lines). Note that all points on a line perpendicular to the dashed line correspond to equal-dose strategies, and that a line closer to the origin corresponds to a smaller number of doses. Clearly, there is a point (p 1 , p2 ), marked by a solid dot, at which R e = (p 1 , p2 ) = 0.85 < 1 so it provides a more effective strategy than (p * 1c , p * 2c ) (equal number of doses but
The time unit is in days.
Similar studies on the role of mixing in vaccination strategies have been conducted with a focus on Economics, in which non-random mixing is referred to as assortative interactions. For example, Galeotti and Rogers 3 studied vaccination strategies when interactions are positive or negative (i.e., assortative or dis-assortative). In this study, two scenarios are discussed, one centralized and other decentralized. The authors aimed to answer the following questions:
(Q1) What information does the government need to efficiently decide the optimal vaccination policy to eliminate an infection? (Q2) Should vaccination be allocated randomly or concentrated in one group? (Q3) How should the planner's action depend on the structure of interactions?
The model considered by Galeotti and Rogers 3 is a system of ODEs for a network consisting of two groups indexed by g ∈ {A, B}. Denote the proportion of the nonimmunized and infected agents in group g by ρ g (t). Assume that a proportion π g Identifying Optimal Vaccination Strategies via Economic and Epidemiological Modeling 5 of group g has immunity. The ODE system for ρ g (t) has the following form:
where ν denotes the per-contact rate of infection, k denotes the number of neighbors that each individual has, δ is the per capita rate at which individuals in the infected state I transition to the susceptible state S, and θ g denotes the probability that a randomly contacted individual in group g is infected, defined by
( We point out that the main results presented in Galeotti and Rogers 3 are for the special case when heterogeneity between the two groups in the transmission rate (ν), connectivity or activity (k), assortativity (β), and group size is ignored. In this paper, we study a more general model by incorporating heterogeneity in these factors.
Our paper is organized as follows. In Sec. 2, we present the generalized model of (1.10) and derive the effective reproduction number R v . We show that R v = 1 is a threshold for the elimination of disease by considering model equilibria and their stability. Section 2.1 focuses on identification of the optimal vaccine allocation. We discuss the results and compare those from the two frameworks (using assortativity or preferential mixing for interactions between groups).
To incorporate heterogeneity, whose importance for disease transmission and control Feng et al. 1 demonstrated in model (1.10) and (1.11), we allow group-specific transmission rate (ν g ), connectivity/activity (k g ), assortativity (β g ), and population size (N g ), g = A, B. Because of the added heterogeneity, model analysis and
Zhao & Feng the identification of optimal immunization strategies is more challenging. But it is important to compare results from simpler (1.10) and (1.11) and more general models, and to determine if the more general model provides any new insights in answering the questions Q1-Q3 mentioned in Sec. 1. As the epidemiologic model includes heterogeneity, this generalization will also permit us to compare these parallel literatures.
In models (1.10) and (1.11), it is assumed that the two groups have identical transmission rates ν, activities k, and assortativities β. If we allow these parameters to differ between groups, the model can be written as
where
which represents the intensity of infection for group g and can differ between groups if k g or ν g are different. All parameters and variables are listed in Table 1.
The pathogen can be eliminated if the disease-free equilibrium, denoted by E 0 = (0, 0), is stable. Let J(E 0 ) denote the Jacobian matrix of (2.1) at E 0 . Then,
where
and Defining
Then H is a non-negative matrix (all elements are non-negative). From the Perron-Frobenius theorem, all eigenvalues of J(E 0 ) have negative real parts if and only if the dominant eigenvalue of H is less than 1. Let R e = r(H) denote the larger eigenvalue of H.
Then
where h ij is defined in (2.6) as functions of model parameters including vaccine coverage π g . Note that R e < 1 if and only if the disease-free steady-state E 0 is locally asymptotically stable; and thus, the pathogen will die out. Using the approach of Reluga et al., 8 we can show that R e < 1 implies global stability of E 0 .
In addition, because all entries of H decrease with π A and π B , we know from Perron-Frobenius Theorem that r(H) (i.e., R e ) will be reduced when π A or π B is increased. Therefore, we can determine the optimal immunization strategy by identifying the combination of (π * A , π * B ) for which R e (π * A , π * B ) = 1 with the smallest number of vaccine doses.
Let N g denote the population size of group g. Then the total number of immunized people, a function of π g , is
(2.8)
Thus, an optimal vaccine allocation (π * A , π * B ) is a solution to the following optimization problem:
Let ∇F = ( ∂F ∂x1 , ∂F ∂x2 ) denote the gradient vector of a scalar function F (x 1 , x 2 ). Note that ∇M = (N A , N B ). An optimal solution (π * A , π * B ) can be obtained by solving jointly the equations
where ω is a constant (Lagrange multiplier). Graphically, a solution to the equations in (2.10) is an intersection of the constraint curve R e (π A , π B ) = 1 (represented by the thick curve) with a contour line given by π
NB (in the unit square), with the smallest c (the blue dashed line). The dot marks the point representing the optimal vaccine allocation (i.e., proportions π * A and π * B of groups A and B immunized to attain R e = 1 using the least amount of vaccine).
We observe in Fig. 2 that, depending on the relative strength of transmission in the two groups (λ A and λ B ) and values for assortativity (β A and β B ), the optimal immunization strategy (π * A , π * B ), marked by the dot, can appear either on the boundary of the unit square (π A , π B ) ∈ [0, 1] 2 or as an interior point with 0 < π * g < 1. For example, in (a) and (b), we have π * A = 0 and π * B = 0, respectively. In (c) and (d), we have π * A = 1 and π * B = 1, respectively. These four cases cover various combinations of β A and β B including both positively and negatively assortative interactions. In (e) and (f), both β A and β B are greater than 1/2, and 0 < π * g < 1 for g = A, B. Other parameter values are:
From λ g = ν g k g /d, the values of ν g , g = A, B, can be determined.
Remark. Unlike the homogeneous case, β A = β B , studied by Galeotti and Rogers, 3 there is no longer a single threshold β = 1/2 to separate optimal immunization strategies for positively or negatively assortative interactions. Also, the optimal allocation for the case β > 1/2 need not be π * A = π * B . Thus, heterogeneity substantially complicates the answer to question Q2.
In Fig. 2, the population sizes for the two groups are equal (N A = N B = 500). Figure 3 illustrates the role of heterogeneous sub-population sizes (N A = N B ) using the same values for other parameters as in Fig. 2(e). Because λ A < λ B and β A < β B , the vaccine coverage is lower in group A than group B. However, when N B > N A , the optimal solution has a higher coverage in group A.
Optimal immunization strategies based on solving the optimization problem (2.9), represented graphically in Fig. 2, identify combinations of immunization proportions
A , π * B ) that reduce R v to 1 with minimal vaccine doses. Additional insights can be obtained by examining the numerical solutions of system (2.1). But note first that the effect of immunization coverage π g on disease prevalence cannot be assessed using the quantity ρ g , which is the fraction infected among unimmunized people. Consequently, even if ρ g is large, the number of infected (I g ) could be very small when if the proportion unimmunized was small. The prevalence can be represented by
Figure 4 illustrates numerical solutions of system (2.1) to compare the effect of various immunization programs on disease prevalence. Other parameter values are chosen to correspond with those in Fig. 2(e). In this case, the optimal strategy to eliminate the pathogen corresponds to the threshold immunization proportions π * A = 0.37 and π * B = 0.71. That is, for all
For the top panel A, we fix π B = π * B = 0.71 and vary π A . We change the plots in the fractions infected in group A (I A /N A , dot-dashed), group B (I B /N B , dotted), and their average (solid) over time. The plots in this panel show that, for π A < π * A , prevalence stabilizes at a positive level, which decreases with π A (see (A1-A3)), and the pathogen is eliminated when π A exceeds π * A = 0.37 (see A4). The plots (B1-B3) show solutions for the same values of R e as those in (A1-A3), respectively, but with different combinations of π g . We observe that, while the average prevalence is similar between panels A and B, the endemic levels in the two groups differ considerably.
We include in this section immunization strategies deduced from the simpler model (1.10) and (1.11) using the same analysis; that is, solution of problem (2.9). First, because groups A and B have identical parameter values except π g , we can simplify the equation R e (π A , π B ) = 1, which from (2.7) yields:
(2.11)
By squaring both sides of the above equation and simplifying further, we obtain the following equation (note that this may introduce an additional solution):
1
which, in the original notation, is
(2.13)
Then, replacing R(π A , π B ) = 1 by G(π A , π B ) = 1 in (2.10), where G is defined in (2.13), we can determine the optimal solution by solving the system:
(2.16)
Using the first equation in (2.15) and eliminating ω, we arrive at the following system, equivalent to (2.15):
(2.17)
Because λ > 0, the first equation in (2.17) implies that for all β = 1/2,
Thus, all optimal solutions satisfying 0 < π * < 1 can be determined by solving the following equation for x:
which can be written as
where
(2.21)
It can be verified that a 2 1 -4a 2 a 0 = 4λ 2 (1β) 2 ≥ 0. Thus, the quadratic equation (2.20) has two real solutions given by
After simplification, we obtain
Recall that x ± = π * A = π * B given in (2.23) are possible solutions only for β = 1/2. In addition, we need to ensure that x ± ∈ [0, 1].
It is easy to show that, when π A = π B = 0, R = λ. This implies that, if λ < 1, E 0 is stable, so there is no need to vaccinate. Thus, we assume that λ > 1. Then, It can also be verified that
Therefore, for 1 2 (1 + 1 λ ) < β < 1, both x -and x + are possible solutions. We then need to determine which is smaller. Note that
That is, x + < x -. It follows that x + is the optimal solution. Figure 5 illustrates the solutions under various conditions. The scenarios presented are based on values of β separated by the following thresholds: . Recall that we assumed that π A ≤ π B . In (c),
x -is the only solution. In (d), both x -and x + are solutions (thus there are two curves corresponding to G = 1), but the optimal solution is π * A = π * B = x + .
The previous sections covered the optimal solution for a central government wishing to vaccinate. However, in the absence of such an authority, people will make their own individually optimal choices. This would almost certainly result in societally suboptimal levels of vaccination, because people do not internalize the negative externalities that they impose on others. Following Galeotti and Rogers, 3 people choose whether or not to be vaccinated before interacting with others. Individuals derive a flow utility of one when healthy and zero when sick. Assuming that people do not have preferences about when they are sick, the utility for each person can be measured by the proportion of time that they are sick. Thus, with cost c > 0 of being immunized, individuals who choose to be vaccinated will ensure themselves a utility of 1c g becauseassuming vaccination is 100% efficacious -they never get sick, but must pay the cost. Alternatively, unvaccinated individuals in group g who spend ρ g proportion of their time sick will have a utility of 1ρ g .
We show in this section that, depending on parameter values, there can be a unique stable interior equilibrium with 0 < π g < 1, g = A, B, or boundary equilibrium points including π A ∈ (0, 1), π B = 1 and π A = 1, π B ∈ (0, 1). For vaccination to be necessary in group g, assume that
) denote an interior equilibrium (i.e., 0 < π * g < 1 for g = A, B). Similar to Galeotti and Rogers, 3
Thus, for π * g to be an interior game equilibrium, the endemic equilibrium (ρ * A , ρ * B ) satisfies For general values of π g , a disease equilibrium ρ * g ∈ (0, 1) of the extended system (2.1) is obtained by setting the right-hand sides of the equations equal to zero, which is equivalent to a positive solution of the following system:
To find an interior game equilibrium π * g ∈ (0, 1), we can solve a similar set of equations as in (3.4) but with ρ * g = c g for g = A, B, and with F (π A , π B ) and G(π A , π B ) being functions of π g :
For ease of presentation, introduce the notation
and let
Lemma 3.1. Let H g be as defined in (3.7). Assume that λg = λ g (1c g ) > 1 and β A + β B = 1. Then a unique interior game equilibrium (π * A , π * B ) exists if the following conditions hold:
(ii)
Moreover, the game equilibrium is given by:
) is a solution of system (3.5), which is linear in 1π A and 1π B . The condition β A + β B = 1 implies that the matrix 2nd Reading August 8, 2019 9:57 WSPC/S0218-3390 129-JBS 1940001
is non-singular. Thus, system (3.5) has a unique solution, 1π * A and 1π * B , from which we obtain:
(3.12)
It can be checked that π * g < 1 if and only if H g /(β A + β B -1) > 0 for g = A, B, which is true if and only if (3.8) or (3.9) holds. It can also be verified that π * g > 0 if and only if H g /(β A + β B -1) < 1, which is true if and only if (3.10) holds.
This completes the proof.
Note that conditions (3.8)-(3.10) only define implicitly regions in the (β A , β B ) plane where the game equilibrium Π * exists. These conditions can be described by the following four lines:
where Q i (1 ≤ i ≤ 3) denote the following constants:
More specifically, the region where Π * exists is the intersection of the regions above (below) the four lines if
To get a better idea what these regions look like, consider the case when the two groups have the same parameter values except c g and β g . In this case, it can be verified that, when β A + β B -1 > 0, the following relations hold: When β A + β B -1 > 0, the directions of inequalities involving Q i are reversed (Fig. 6).
We observed in Fig. 6 that Π * does not exist for all β g , and that the size of the region for existence increases as c A and c B become similar. The parameter values used are the same as before except that ν A = ν B = 0.12, k A = k B = 4. In (a) c A = 0.12 > c B = 0.1, and in (b) c A = 0.1 < c B = 0.13.
Remark. In the special case, when all parameters are identical in the two groups (except π g ), as in the model of Galeotti and Rogers, 3 Q i = 1 for i = 1, 2, 3. Thus, all four lines defined in (3.13) are the same as β A + β B = 1 (i.e., 2β = 1 when β A = β B ). Therefore, Π * exists for all parameter values except β = 1/2.
For the stability of game equilibrium Π * , we can examine the derivative of the best response of group B with respect to group A's immunization rate, as done by Galeotti and Rogers, 3 and identify conditions under which
(3.15) Lemma 3.2. A necessary condition for the interior game equilibrium Π * = (π * A , π * B ) to be stable is the following inequality: Proof. For ease of notation, let ρ A (π A ) and dρ A (π A )/dπ A denote the following quantities:
(3.17)
Differentiating the equation G(ρ A , ρ B ) = 0 in (3.5), evaluating at ρ * A = ρ A (π * A ), and noting that ρ B (π A , π * B (π A )) = c B when π A is sufficiently close to π * A , we have
2), we have
(3.20)
To use the condition P B > -1, we show first that the denominator in (3.20) is positive. It suffices to note that
and thus,
From (3.20), using P B > -1 and the fact that the denominator of (3.20) is positive, we know that a necessary condition for Π * to be stable is the following inequality: Remark. The necessary condition (3.16) for the stability of the interior game equilibrium is much simpler in the model considered in Galeotti and Rogers. 3 This is because, when λ g = λ, c g = c, β g = β, (3.16) simplifies to:
which can be true if and only if 1 -2β < 0 (i.e., β > 1/2). Thus, a stable interior game equilibrium can be stable only for positive assortative interactions. However, this may not be the case in the extended model (2.1).
Although inequality (3.16) does not provide a simple condition in terms of assortativity, we can examine it numerically (Fig. 7). The contour curve of the function L(β A , β B ), which has the property that L > 0 is equivalent to the inequality (3.16), where
When the c g change, the game equilibrium Π * may change. Figure 8 shows the dependence of π * A and π * B on costs c g . All other parameter values are the same as in Fig. 2(e).
We observe from Figs. 8(a) and 8(b) that, for given c B , π * A decreases with c A (this is intuitive inasmuch as a higher cost would lead to lower proportion immunized) while π * B increases with c A . Symmetric properties can also be observed (i.e.,
known to influence disease transmission and control via immunization, and analyze the resulting model using an approach devised for epidemiological models. Our model (2.1) is capable of generating heterogeneous interior optimal vaccination strategies, and can produce qualitative and quantitative outcomes very different from the simpler model (1.10). For example, Galeotti and Rogers 3 pointed out that an interior equilibrium solution (π * A , π * B ) of model (1.10) can only happen in symmetric form; i.e., π * A = π * B ∈ (0, 1), and that any asymmetric equilibrium will occur on a boundary; i.e., either π A = 0 or π B = 1 (assuming that π A < π B ). In our extended model (2.1), however, there are asymmetric interior equilibria (see Figs. 2(e) and 2(f)).
We show that optimal strategies in the new model are no longer determined solely by the nature (positive or negative) of assortative interactions. Consequently, the answers to questions Q1-Q3 are complicated by additional threshold conditions. Both analytic and numerical analyses of the new model are conducted to investigate how vaccination strategies may depend on both disease transmission parameters and costs associated with vaccination and illness. Optimal vaccination strategies under both centralized vaccination programs and decentralized games are discussed.
An interesting question is how are the models formulated by economical modelers (such as the model (1.10) in Galeotti and Rogers 3 ) and the model formulated by epidemiological modelers (such as model the simplest SIS or SIR model) related? If they are significantly different, what are the underlying assumptions in addition to those already pointed out? In Appendix A.1, we tried to recover from the equations for proportions (ρ g ) in (1.10) the equations in terms of the numbers of susceptible (S g ), infected (I g ), and immune (R g ). The result seems to suggest that some biologically unreasonable assumptions have been made in model (1.10) (see the remarks at the end of Appendix A.1). One of the key assumptions is that people do not gain immunity from infections while the immunity induced from vaccination does not wane. However, for almost all vaccine preventable infectious diseases, if vaccination induces an immune response, so must infection. If we allow vaccine-induced and naturally-acquired immunity differ (e.g., in the degree or duration of protection), there should be two separate classes, e.g., V (for vaccinated and immune) and R (for recovered with temporary immunity). This is demonstrated in Appendix A.2. Clearly the modified model (A.8) based on a SIR model is more biologically reasonable than model (2.1) based on a SIS model. One of the main differences between the two systems will be in the analysis that involves the endemic equilibrium, such as the effect of vaccination on disease prevalence or the game equilibrium for the decentralized game. More detailed analyses for system (A.8) or a more general model will be presented in future work.
Analysis of System (2.1) suggests that the optimal vaccination strategies (π * A , π * B ) under centralized and decentralized cases may depend on the cost function. If so, this will illustrate the importance of melding Economics and Epidemiology. We have begun studying this and will publish results elsewhere.
Zhao & Feng
August 8, 2019 9:57 WSPC/S0218-3390 129-JBS