key: cord-0214897-29dydv2e
authors: Jhun, Bukyoung; Choi, Hoyun
title: Efficient vaccination strategies to reduce the epidemic mortality in population with heterogeneous fatality rate
date: 2021-09-08
journal: nan
DOI: nan
sha: bd929adecaaae69e0622e044117fa4bbbcfe61bb
doc_id: 214897
cord_uid: 29dydv2e

An insufficient supply of effective SARS-CoV-2 vaccine in most countries demands an effective vaccination strategy to minimize the damage caused by the disease. Currently, many countries vaccinate their population in descending order of age to minimize the deaths caused by the disease; however, the effectiveness of this strategy needs to be quantitatively assessed. We employ the susceptible-infected-recovered-dead (SIRD) model to investigate various vaccination strategies. We constructed a metapopulation model constructed from empirical human contact and fatality of SARS-CoV-2 and investigated vaccination strategies. We found that the age-based strategy, which is currently employed in many countries, is more effective when the basic reproduction number is high and vaccination supply is low, but the rate-based method outperforms the age-based strategy when there is sufficiently high supply of vaccine. Simulated annealing is performed to find the optimal vaccination strategy. We identified a first-order phase transition in the vaccination strategies which is characterized by a discontinuous transition in the optimal strategy and the hysteresis. This phase transition implies that combining the age-based and rate-based strategy is ineffective in reducing the number of deaths. These conclusions are valid even when the heterogeneous degree distribution of human contact is considered.

The spreading process in complex systems such as networks [1] [2] [3] [4] [5] and metapopulation [6] [7] [8] [9] [10] has been an active field of research for modelling many physical and social phenomena [11] [12] [13] . These have included opinion formation in social groups [13] [14] [15] [16] , the spread of epidemic diseases [6, [17] [18] [19] [20] [21] , and innovations [22, 23] . A plethora of data [24] [25] [26] on human mobility, collaboration, contagion of epidemic disease, and temporal contacts, all of which were previously unavailable to researchers, now enable effective research into various dynamic processes in social systems. Extensive research devoted to the spreading processes has provided quantitative analyses for policy-making and especially, in the public health domain. Moreover, the spreading process has provided a deeper understanding of phase transitions and critical behaviors, such as the effect of structural heterogeneity on epidemic thresholds [6, 19, 27] and the hybrid phase transition induced by cascades [28] [29] [30] [31] .

One of the most important topics in mathematical epidemiology is vaccination strategy, which has been extensively studied with various epidemic models [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] . If an individual is vaccinated, that individual will not turn into the infected state even if it is in contact with an infected individual. Moreover, when a sufficient fraction of individuals in a system are vaccinated, the infection is unable to spread throughout the system, and the epidemic state is eliminated by the vaccination. This effect is called herd immunity. Vaccination strategies frequently aim to achieve herd immunity with the smallest number of vaccine shots.

The SARS-CoV-2 pandemic is ongoing worldwide, and has caused more than four million deaths. Due to the development of effective vaccines for the disease, the elimination of the disease can be anticipated. However, in most countries, and especially in the developing countries, the number of vaccine shots * jhunbk@snu.ac.kr available is less than the total population [47] . Therefore, it is important to formulate a vaccination strategy that minimizes the damage caused by the disease, such as the number of deaths, with the limited supply of vaccine available. Currently, many countries are vaccinating their populations in descending order of age, since the infection fatality rate (IFR) for the COVID-19 increases with age [48] [49] [50] [51] [52] [53] [54] . However, the effectiveness of this strategy needs to be quantitatively assessed.

Here, we employ the susceptible-infected-recovered-dead (SIRD) model, which is a minimal model to study the epidemic mortality. We evaluate the effectiveness of age-based and rate-based vaccination strategies in a metapopulation model with empirical data obtained from previous research (see Sec. III A for details). We find that the age-based strategy is more effective than the rate-based strategy for a high basic reproduction number and low vaccination supply, but the rate-based strategy outperforms the age-based strategy when sufficiently large amount of vaccine is supplied. Simulated annealing is implemented to find the globally optimal vaccination strategy. We find that there is a first-order phase transition in the vaccination strategy characterized by an abrupt transition in the average age of the vaccinated population and hysteresis. Further, the calculation is performed in a metapopulation with a heterogeneous degree distribution to find that the conclusions derived from a metapopulation model with a homogeneous degree distribution are valid when a heterogeneity of the degree distribution is considered. This paper is organized as follows: First, we introduce the SIRD model in Sec. II. Next, in Sec. III A, we explain the empirical data used to construct the metapopulation model for SARS-CoV-2 epidemics. In Sec. III B, we investigate the efficiency of vaccination strategies in the metapopulation model constructed from the empirical data, and in Sec. III C, we study the phase transition of the vaccination strategy. In Sec. III D, we show that the conclusions from the previous sections are valid even when the heterogeneous degree distribution of human contacts is considered. A summary and final remarks are presented in Sec. IV.

The susceptible-infected-recovered (SIR) model is a minimal model of epidemic spreading and the most extensively studied model both in complex networks [17, 18, 27, 55] and in the metapopulation model [6] [7] [8] [9] [10] , together with its variants [28, [56] [57] [58] [59] [60] [61] . In the SIR model, each individual is in either the susceptible (S), infected (I), or recovered (R) state. A susceptible individual can turn into the infected state if it comes into contact with an infected individual. The rate of transition is proportional to the number of infected individuals in contact. Infected individuals eventually turn into the recovered state at a constant rate. A recovered individual obtains immunity and does not turn into the infected state again. In the real world, actual infectious diseases exhibit more complex stages; however, they can be reduced roughly to the SIR model. Therefore, conclusions from a study employing the SIR model are applicable to actual infectious diseases.

Vaccination strategy is one of core topics in mathematical epidemiology; therefore, considerable research has been devoted to the subject [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] . The objective of vaccination strategies in the SIR model is to minimize the total number of individuals affected by the disease, which can be measured by the number of recovered individuals when the infection vanishes, with limited vaccination resources. However, one of the most important objectives of vaccinations in the real world is to minimize the total number of deaths caused by a disease. Since recovery and death are not distinguished in the SIR model, the model cannot be used to study vaccination strategies relating to such a purpose. At this point, we employ the SIRD model, which is a minimal model that distinguishes recovery and mortality [62] [63] [64] [65] . By the same spirit of the SIR model, it should be possible to reduce the complex stages exhibited by the actual disease to the SIRD model; hence, the results of this research should be valid for actual infectious diseases.

In the SIRD model, similar to the SIR model, each individual is in either a susceptible (S), infected (I), recovered (R), or dead (D) state. If a susceptible individual and an infected individual are in contact, the susceptible individual is turned into the I state at rate η (it turns with probability η∆t in an infinitesimal time step ∆t). If the S individual is in contact with n infected individuals, the rate becomes nη. Any individual from subpopulation α (such as age group) that is in the I state turns into the R state at rate µ, or into the D state at rate κ α µ. The rate equation for the SIRD model is, therefore,

where η is the contagion rate, µ is the recovery rate, and κ α is a non-negative parameter associated with subpopulation α that controls the IFR of the disease. The IFR is defined as the ratio of deaths from a certain disease compared to the total number of people infected with the disease. If an individual from subpopulation α is infected, the individual turns into R state or D state with probability ratio κ α : 1. Therefore, the IFR of the subpopulation α is

which is always between zero and one for non-negative κ α . If the contagion rate of an epidemic disease is higher than a threshold, an epidemic outbreak that starts from an infinitesimal fraction of individuals can spread over an entire system. One way to parametrize epidemics is through the basic reproduction number [66] [67] [68] , which is usually denoted as R 0 . The basic reproduction number is the anticipated number of contagions caused by a single infected individual until the individual becomes unable to spread the infection anymore, due to either recovery or death. If R 0 > 1 at t = 0, on average, a single infected individual will spread the infection to more than one individual; thus, the infection will spread. Therefore, the threshold of an epidemic can be expressed as R 0 = 1. Initially, only an infinitesimal proportion of the individuals are infected; therefore, the basic reproduction number can be expressed as [69, 70] 

whereΛ is the largest eigenvalue of the modified contact matrix M αβ /(1 + κ β ) of the metapopulation. M αβ is the contact matrix (see Sec. III B for details) of the metapopulation model, and the term (1+κ β ) accounts for the fact that in subpopulation β, the rate in which infected state turns into any other state is µ(1 + κ β ) instead of µ.

The contact matrix M αβ is defined as the average contacts that an individual in group α has with the individuals in group β. Contact data between each age group in United States has been studied by means of survey [71] . The contact matrix of the United States population for each age group is illustrated in Fig. 1(a) . The population is divided into 9 groups: aged 0-9, 10-19, ..., 70-79, and above 80. Contact within the similar age group is disproportionately intense, and contacts between teenagers exhibit the highest strength.

The IFR of SARS-CoV-2 is an active topic of research [48] [49] [50] [51] [52] 54] . One challenge in the estimation of the fatality rate of an epidemic disease is that if there are large number of unconfirmed infections, the fatality rate can be overestimated. The IFR of SARS-CoV-2 in the US cannot be accurately estimated from the case fatality rate (CFR), which is the ratio of deaths caused by a disease to the total number of people diagnosed with the disease, because of the large disparity between the number of people infected by the disease and the number of people diagnosed with the disease. The number of people infected by the disease can be estimated by the seroprevalence [72] [73] [74] , which is the percentage of people in a population who have antibodies of a certain disease that suggests their exposure to the infectious agent. The seroprevalence of SARS-CoV-2 in the US varies significantly from county to county and is up to 44% in some regions while the rate of confirmed cases in the region is not [75] [76] [77] . To this end, instead of using the fatality data from the US, we use employ up-to-date estimation of the IFR in South Korea [53] . Intensive contact tracing and testing, often referred to as the Kquarantine model [78] [79] [80] [81] [82] [83] [84] , together with the small number of cases, results in an accurate estimation of the IFR. The estimated IFR for each age group in South Korea is illustrated in Fig. 1(b) . This composition of population results in a high epidemic mortality rate: if the population consists solely of young people, almost no people will die because of their extremely low IFR, and if there are only senior individuals in the population, the epidemics cannot spread through the population because of their low contact rate.

A metapopulation model consists of interacting subpopulations, which are often, but not necessarily, spatially structured. The subpopulations are assumed to be well-mixed. For epidemic studies in the metapopulation model, we track the density of epidemic states in each subpopulation instead of tracking epidemic states of each individual. The density of states evolves due to the interactions among subpopulations, including interactions within the same subpopulation. Because metapopulation models have lower dimensions compared to networks, it allows more exhaustive studies on the spread of epidemic diseases. The epidemic equation for the SIRD model in the metapopulation model is

where P S α , P I α , P R α , and P D α are the probability that an individual in group α is in the S, I, R, and D state, respectively.

Initially, an infinitesimal fraction, n 0 = 10 −8 of each group α of the population, is in the I state, and all the rest of the population 1 − n 0 is in the S state. The differential equations are then solved by fourth order Runge-Kutta method [85, 86] until the total proportion of infected individuals α β M αβ P I β becomes less than a certain threshold 10 −12 . We then calculate the total proportion of the dead population α β M αβ P D β . We implement four vaccination strategies: age-based, ratebased, simulated annealing, and random vaccination strategies. Since the IFR of the SARS-CoV-2 monotonically increases with age, in the age-based method, we vaccinate the population in descending order of age ( Fig. 1(c) ). In the ratebased method, infinitesimal amount of vaccine is iteratively given to the age group with the highest contact rate with unvaccinated individuals until the total amount of vaccine is distributed. The contact rate of age group α with unvaccinated individuals is

where v β is the fraction of vaccinated individuals in the group β, and the value is recalculated at each iteration. This differs from the contact rate of age group α with any individual in the population which is β M αβ . The age group vaccinated by the rate-based strategy is illustrated in Fig. 1(d) . At first, the teenagers who have the highest contact rates are primarily vaccinated. As the vaccination progresses, because teenagers mostly interact with each other (Fig. 1(a) ), large fraction of the individuals they contact with is already vaccinated. As a result, the contact rate of teenagers with unvaccinated individuals rapidly decreases, and the population aged 20-29 are then vaccinated. The random vaccination strategy homogeneously vaccinates the population regardless of their age. When the fraction of the population affected by the epidemic, which is the fraction of the R and D state at the end of the dynamic, is less than 10 −4 , we assume that herd immunity is achieved, and the epidemic is eliminated by vaccination. The point at which herd immunity is achieved is depicted as a dot.

(d-f) The fraction of the population in each group vaccinated by the optimal strategy is found by simulated annealing. The dashed white line is the threshold for herd immunity. In (d), the method involves vaccinating the teenagers, and shows no complex behavior. However, there is an abrupt transition of strategy in (e) and (f). For a low vaccine supply, the method involves vaccinating the senior population, but for sufficiently large vaccine supply, the young population is primarily targeted for vaccination. This abrupt transition of the strategy indicates a first-order phase transition.

The simulated annealing is a probabilistic optimization algorithm inspired by spin glass [87] . We implemented a modified version of the simulated annealing technique to find the globally optimal vaccination strategy. First, we start with random vaccination strategy with the given amount of vaccine supply. We set this strategy as a provisional solution. Then we calculate the mortality rate of a trial strategy, which is perturbed from the provisional solution by a small amount, while keeping the vaccine supply of the total population constant. If the mortality rate of the trial strategy is smaller than that of the provisional solution, we replace the provisional solution with the trial strategy. Otherwise, we replace the provisional solution with the trial strategy with probability exp (−1/T ), where T is the temperature of the algorithm. At the beginning, the temperature is set at T = 2. We iterate this process for n iter = 10 6 times, while the temperature is dropped by a factor f iter = 1 − 10 −5 at each step. The resulting provisional solution is the optimal vaccination strategy, given that n iter is sufficiently large and f iter is sufficiently close to one.

We performed numerical calculations using the contact matrix and the IFR data explained in the previous section. The mortality rates resulting from the strategies are illustrated in Fig. 2(a-c) for various basic reproduction numbers R 0 s. The contagion rate η corresponding to each R 0 is calculated by Eq. (5). If the fraction of population in either the R or D state at the end of the epidemic process is less than 10 −4 , we assume that herd immunity is achieved and that the epidemic is eliminated by the vaccination. The age-based strategy is advantageous over the rate-based strategy for large R 0 and a small vaccine supply; however, for sufficiently large vaccine supply, the rate-based strategy outperforms the age-based strategy. For a large vaccine supply, the age-based strategy is even less efficient than the random vaccination strategy. This is because the senior population targeted by the age-based strategy has low contact rate, therefore, vaccinating them is ineffective in obtaining herd immunity. The optimal strategy found by the simulated annealing exhibits mortality rate almost identical to either the age-based or rate-based strategies. This suggests that the globally optimal strategy is close to one of these two strategies. The population targeted by the optimal strategy is illustrated in Fig. 2(d-f) . The dashed white line in the figure marks the point at which herd immunity is achieved by the optimal strategy. There is a sharp transition of the strategy in Fig. 2(e, f) , suggesting the existence of a discontinuous phase transition in the vaccination strategy.

The difference between the mortality rates of age-and rate- The difference between the mortality rate resulting from age-based and rate-based strategies. The dashed line represents the boundary between the regions where the age-and rate-based strategies are advantageous in reducing the mortality rate. The age-based strategy reduces more deaths than the rate-based strategy when the basic reproduction number R 0 is high and the level of vaccine supply is low. However, as the vaccine supply becomes higher, the ratebased strategy outperforms the age-based strategy.

based strategies for various vaccine supply and R 0 is illustrated in Fig. 3 . The contagion rates ηs corresponding to R 0 s are calculated by Eq. (5). When R 0 < R * 0 = 1.12, the ratebased strategy outperforms the age-based strategy for all levels of vaccine supply. When R 0 > R * 0 , the age-based strategy is advantageous for low vaccine supply; however, when the vaccine supply is above a critical value, which is depicted as a dashed line in the figure, the rate-based strategy reduces more deaths than the age-based strategy. The solid line shows where the both strategies achieve the herd immunity hence the difference between the strategies becomes irrelevant.

In this section, we further investigate the nature of phase transition suggested by Fig. 2(e) and (f). To find locally optimal solutions, we use the zero-temperature simulated annealing, which is analogous to the gradient-descent method. Initially, the provisional solution is set to either the age-based or rate-based strategy. Then, we perturb the provisional solution by decreasing the vaccination rate of group α by δv/P (α) and increasing the vaccination rate of group β by δv/P (β), where δv is a small number and P (α) is the fraction of the group α in the population. This way, the total vaccination rate of the entire population remains constant. Among perturbed solutions, if any solution results in a smaller mortality rate, we replace the provisional solution with the perturbed solution that results in the smallest mortality rate. Otherwise (i.e. if all the perturbed solutions result in larger mortality rates than the provisional solution), we have achieved a locally optimal solution, hence we terminate the process.

We searched for a locally optimal solution in the vicinity of the age-based and rate-based strategy in the metapopulation model. We illustrate the average age of the vaccinated population for the locally optimal strategies in Fig. 4 . This value serves as the order parameter of the phase transition. The two results coincide at low levels of vaccine supply, suggesting that there is a single globally optimal strategy in the system. However, for sufficiently large vaccine supply, the two results are distinct, i.e., the system exhibits hysteresis. This implies that there are two locally optimal strategies in the system. The average age of the vaccinated population changes discontinuously at the transition point. The abrupt change of the average age of vaccinated individuals and hysteresis suggest that the phase transition is of first-order. The hysteresis of this phase transition also implies that the moderate strategy of the age-and rate-based strategy is less effective than either strategies. The vaccination rate of the moderate strategy is v mod

and v rate α are vaccination rates of the age group α in the ageand rate-based strategy, respectively. The performance of the moderate strategy for various levels of vaccine supply is depicted in Fig. 5(a) . The mortality rate as a result of the moderate strategy with r = 0.5 is represented by the dashed black line. The moderate strategy is never more effective than both of the age-and rate-based strategies, and in some region, it is less effective than any of the two strategies. Illustrated in Fig. 5(b) is the performance of the moderate strategy versus the ratio r. The basic reproduction number R 0 = 1.3 and the total vaccine supply ranges from 8% to 18%. There is a barrier of mortality rate between the age-based and rate-based strategies, which is the origin of the hysteresis. The mortality rate of the rate-based strategy (r = 0) decreases faster than the age-based strategy (r = 1) and achieves zero mortality rate at a lower level of vaccine supply.

It was reported that the degree distribution of contact between individuals is not homogeneous, but follows a negative binomial distribution NB (r, p) with r ≃ 0.36 [88] . The parameter p β of the age group β can be determined by the average degree k β = α M αβ : p β = 1 − k β / r + k β . Then the differential equations (Eqs. (6)-(8)) are replaced by

where P (α) d (k) is the degree distribution of the age group α. P S αk , P I αk , P R αk , and P D αk are the probability that an individual with degree k in group α is in the S, I, R, and D state, respectively. We track P S αk , P I αk , P R αk , and P D αk for each α and k = 0, 1, · · · , 99. Initially, an infinitesimal fraction n 0 = 10 −8 5. (a) The mortality rate as a result of age-based, rate-based, and the moderate strategies. The moderate strategy is a mixed strategy of the age-and rate-based strategies with r = 0.5. The moderate strategy is never an optimal strategy and less efficient than both of the strategies in some interval. (b) The mortality rate of the moderate strategy of age-and rate-based strategy for R 0 = 1.3 and vaccine supply, from top to bottom, 8%, 10%, 12%, 14%, 16%, and 18%. There is a barrier of mortality rate between the age-and rate-based strategy, which is the cause of the hysteresis of the phase transition. The mortality rate of the rate-based strategy drops faster than the rate-based strategy as the level of vaccine supply increases, causing the crossover between the two strategies.

of each group (α, k) of the population is in the I state and all the rest of the population 1 − n 0 is in the S state. Then, the differential equations are solved until the total fraction of the infected individuals α β k M αβ P (β) d (k) P I βk becomes less than the threshold 10 −12 .

The mortality rate resulting from the random, age-based, and rate-based strategies in the metapopulation model with heterogeneous degree distribution is illustrated in Fig. 6 . The conclusions from the previous sections hold even when the heterogeneity of the degree distribution is considered. A large R 0 gives the advantage to the age-based strategy relative to the rate-based strategy; however, the rate-based method outperforms the age-based method for sufficiently large vaccine supply. The rate-based strategy achieves herd immunity at a lower vaccination rate. The simulated annealing was not performed due to the excessively large computational cost required for the heterogeneous degree distribution.

In summary, we employed the SIRD model to investigate the effectiveness of vaccination strategies to minimize the mortality rate in a metapopulation. We constructed the metapopulation model from a survey of human contacts and an up-to-date estimation of the IFR for SARS-CoV-2. We investigated the age-based, rate-based, and simulated annealing strategies, together with random vaccination as a benchmark. The age-based strategy vaccinates the population in descending order of age, and the rate-based strategy vaccinates first the age groups with the highest contact rate. As a result, the rate-based strategy primarily vaccinates the population of ages 10-29 because of their high contact rate. Simulated annealing was employed to find globally optimal vaccination strategies. The age-based strategy is more effective than the rate-based strategy for a high basic reproduction number with a low level of vaccine supply; however, the rate-based strategy outperforms the age-based strategy when a high level of vaccine supply is available and eventually achieves herd immunity at a lower vaccine supply. We found that there is a first-order phase transition in the optimal vaccination strategy which is characterized by a discontinuous transition of globally optimal strategy and hysteresis. The hysteresis of this phase transition implies that combining these two strategies is ineffective in reducing the mortality rate of the epidemic disease. Therefore, if the age-based or rate-based strategy has been implemented, changing the strategy to the other is likely to be ineffective. These conclusions are shown to be valid when the heterogeneous degree of distribution of human contact is considered.

In this study, the vaccine is distributed before the epidemics The rate of recovery is set µ = 1. If the fraction of the population affected by the infection, which is the density of individuals that are in either the R or D state at the end of the epidemics, is below 10 −4 , we assume that the epidemic is eliminated by the vaccination and herd immunity is achieved. Such a point is depicted as a dot in each plot. For large R 0 , the age-based strategy outperforms the rate-based strategy for a low level of vaccine supply; however, for a sufficiently high level of vaccine supply, the rate-based method outperforms the age-based strategy and achieves the herd immunity at lower vaccination rate.

occur with an infinitesimal fraction of infected individuals. These are simplified assumptions because, in the real-world, a nonvanishing fraction of individuals are infected by the time vaccine is deployed, and the vaccine cannot be administered instantly. Nonetheless, epidemics started from a finite fraction of infected individuals often exhibit similar qualitative behaviors with the ones started from infinitesimal infectious seeds [17, 28, [89] [90] [91] [92] . The efficiency of each vaccination strategy will be changed if the vaccination is temporally deployed in parallel with the epidemics; however, the qualitative results of this study will still be valid because the barrier of mortality rate in the efficiency landscape will not disappear when such perturbation is present.

In conclusion, the effectiveness of vaccination strategies is closely related to the amount of vaccine available. Hence, the quantity of vaccine supply should be estimated before the design of the vaccination strategies. Precise estimation of the contact matrix, basic reproduction number, and the IFR of the population is also important. In the survey data used in this paper, all types of contacts are treated equally. However, the contagion rate of disease between individuals who live in the same house, work in the same place, or shop in the same grocery store should differ from each other. If more accurate contagion tree data of the disease is collected and implemented, the relative strength of such interactions will be taken into account. Although the effectiveness of the strategies at specific vaccination rates will be modified if the precision of the dataset is improved, because the nature of the phase transition is invariant under perturbations, the phenomena observed in this paper and the conclusions should still be valid.

Proc. Twelfth ACM Int. Conf. Web Search Data Min

2020 IEEE Intl Conf Parallel Distrib

Further Notes on the Basic Reproduction Number

Readings Comput. Vis

2013 Winter Simulations Conf

Proc. 5th ACM Conf

This research was supported by the NRF, Grant No. NRF-2014R1A3A2069005.