# Abrupt transition of the efficient vaccination strategy in a population with
heterogeneous fatality rates

Bukyoung Jhun<EMAIL_ADDRESS>CCSS, CTP and Department of Physics and
Astronomy, Seoul National University, Seoul 08826, Korea  Hoyun Choi  CCSS,
CTP and Department of Physics and Astronomy, Seoul National University, Seoul
08826, Korea

###### Abstract

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 (i.e. descending order of fatality rate) 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 with heterogeneous contact and fatality rates and
investigated the effectiveness of vaccination strategies to reduce epidemic
mortality. We found that the fatality-based strategy, which is currently
employed in many countries, is more effective when the contagion rate is high
and vaccine supply is low, but the contact-based method outperforms the
fatality-based strategy when there is a sufficiently high supply of the
vaccine. We identified a discontinuous transition of the optimal vaccination
strategy and path-dependency analogous to hysteresis. This transition and
path-dependency imply that combining the fatality-based and contact-based
strategies is ineffective in reducing the number of deaths. Furthermore, we
demonstrate that such phenomena occur in real-world epidemic diseases, such as
tuberculosis and COVID-19. We also show that the conclusions of this research
are valid even when the complex epidemic stages, efficacy of the vaccine, and
reinfection are considered.

> The effectiveness of vaccination depends highly on the choice of the
> individuals to vaccinate, even if the same number of individuals are chosen.
> Therefore, effective vaccination strategies have been a central topic of
> research in mathematical epidemiology and provided quantitative analysis to
> inform policy-making in the public health domain. In this study, we
> investigated the effectiveness of vaccination strategies to reduce epidemic
> mortality in a population where the fatality rate varies among groups of
> individuals by constructing a metapopulation model with heterogeneous
> contact and fatality rates. We show that the effectiveness of vaccination
> strategies is closely related to the amount of vaccine available. When the
> vaccine supply is low, vaccinating individuals in the descending order of
> fatality rates is effective; while when the vaccine supply is sufficiently
> high, vaccinating individuals with high contact rates outperforms the
> fatality-based strategy. By employing simulated annealing, we identify an
> abrupt transition of the optimal strategy and path-dependency analogous to
> hysteresis in statistical physics. This transition suggests that combining
> the fatality- and contact-based strategies is less effective than either
> strategy; therefore, if a country has been employing a specific vaccination
> strategy, it is inadvisable to convert from that strategy to the other. We
> also show that such phenomena occur in the vaccination of real-world
> epidemic diseases, such as tuberculosis and COVID-19.

##  I Introduction

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
modeling many physical and social phenomena [11, 12, 13]. This research has
included opinion formation in social groups [13, 14, 15, 16], the spread of
epidemic diseases [17, 18, 6, 19, 20, 21], and the diffusion of innovations
[22, 23, 24, 25]. Current access to a plethora of data [26, 27, 28] on human
mobility, collaboration, the 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, especially in the public health domain. Moreover, study of
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, 29] and the hybrid phase transition induced by
cascades [30, 31, 32].

One of the most important topics in mathematical epidemiology is vaccination
strategy, which has been extensively studied with various epidemic models [33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]. If an individual is
vaccinated for certain epidemic disease, that individual acquires immunity to
the disease. Actual vaccines have less than perfect efficacy, which means that
there is a small probability that a vaccinated individual can be infected by
the disease (i.e., a vaccine breakthrough). It is often modeled that
vaccinated individuals do not turn into the infected state even in contact
with infected individuals. In such a model, 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 five
million deaths to date. Due to the development of effective vaccines for the
disease, the epidemic damage of the disease can be greatly reduced. However,
in most countries, especially developing countries, the number of vaccine
shots available is less than the total population [48]. 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
vaccines 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 [49, 50, 51, 52, 53, 54, 55].
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 epidemic mortality. We evaluate the effectiveness of
fatality-based and contact-based vaccination strategies in a metapopulation
model with heterogeneous contact and fatality rates. We find that the
fatality-based strategy is more effective than the contact-based strategy for
a high contagion rate and low vaccination supply, but the contact-based
strategy outperforms the fatality-based strategy when a sufficiently large
amount of vaccine is available. Simulated annealing is implemented to find the
globally optimal vaccination strategy. We find that there is a discontinuous
transition of the optimal strategy and path-dependency analogous to
hysteresis. Further, we demonstrate that these phenomena occur in the
vaccination of real-world epidemic diseases, such as tuberculosis (TB) and
COVID-19.

This paper is organized as follows: First, we introduce the SIRD model in Sec.
II. Next, in Sec. III.1, we introduce the synthetic metapopulation model
constructed for this research, which has heterogeneous fatality and contact
rates. In Sec. III.2, we study the transition and path-dependency of the
vaccination strategy. In Sec. III.3, we show that such a transition occurs in
real-world epidemic diseases, such as tuberculosis and COVID-19, and in Sec.
III.4, we demonstrate that the results of this research are valid even when
complex details of the epidemic diseases are considered. A summary and final
remarks are presented in Sec. IV.

##  II Susceptible-infected-recovered-dead (SIRD) model

The susceptible-infected-recovered (SIR) model is a minimal model of epidemic
spreading and the most extensively studied model both in complex networks [18,
17, 29, 56] and in the metapopulation model [6, 7, 8, 9, 10], together with
its variants [30, 57, 58, 59, 60, 61, 62, 63, 64, 65]. In the SIR model, each
individual is in either the susceptible (S), infected (I), or recovered (R)
state. A susceptible individual can turn into an infected state if it comes
into contact with an infected individual. If a susceptible individual and an
infected individual are in contact, the susceptible individual is turned into
the I state at rate $\eta$ (it turns with probability $\eta\Delta t$ in an
infinitesimal time step $\Delta t$). If the S individual is in contact with
$n$ infected individuals, the rate becomes $n\eta$. Infected individuals
eventually turn into the recovered state at a constant rate $\mu$. A recovered
individual obtains immunity and does not turn into the infected state again.
In actual epidemic diseases, there is a probability of reinfection whose
effects on the results of this research are discussed in Sec. III.4.

Vaccination strategy is one of the core topics in mathematical epidemiology;
therefore, considerable research has been devoted to the subject [33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]. 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. Because
recovery and death are not distinguished in the SIR model, it cannot be used
to study vaccination strategies related to such a purpose. At this point, we
employ the SIRD model, which is a minimal model that distinguishes recovery
and mortality [66, 67, 68, 69].

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. The contagion
occurs identically as in the SIR model. Any individual from subpopulation
$\alpha$ (such as an age group) that is in the I state turns into the R state
at rate $(1-\kappa_{\alpha})\mu$ or into the D state at rate
$\kappa_{\alpha}\mu$. The rate equation for the SIRD model is, therefore,

$\displaystyle\mathrm{S}+\mathrm{I}\overset{\eta}{\rightarrow}\mathrm{I}+\mathrm{I}\,,$
(1)
$\displaystyle\mathrm{I}\overset{(1-\kappa_{\alpha})\mu}{\rightarrow}\mathrm{R}\,,$
(2)
$\displaystyle\mathrm{I}\overset{\kappa_{\alpha}\mu}{\rightarrow}\mathrm{D}\,,$
(3)

where $\eta$ is the contagion rate, $\mu$ is the recovery rate, and
$\kappa_{\alpha}$ is the IFR of subpopulation $\alpha$. If an individual from
subpopulation $\alpha$ is infected, the individual turns into R state or D
state with probability ratio $(1-\kappa_{\alpha}):\kappa_{\alpha}$. We assumed
that the three processes (contagion, recovery, and death) occur independently
at constant rates. This assumption reasonably describes the pathology of each
individual; however, complex social interventions such as quarantine and
social distancing that depend on the number of epidemic cases and mortality
can complicate the process.

##  III Results

###  III.1 Fatality- and contact-based strategies

Figure 1: (a–b) The mortality rate as a result of various vaccination
strategies. The contagion rate is (a) $\eta=0.05$ and (b) $\eta=0.4$, while
the recovery rate is normalized $\mu=1$. When the fraction of the population
affected by the epidemic (i.e., the fraction of 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. When the contagion rate is low, the
contact-based strategy is more effective regardless of the vaccination rate.
When the contagion rate is high, the fatality-based strategy is more effective
at a low vaccination rate; however, the contact-based strategy outperforms the
fatality-based strategy when the vaccine supply is sufficiently high,
achieving herd immunity at a lower vaccination rate. (c) The difference
between the mortality rates resulting from fatality-based and contact-based
strategies. The fatality-based strategy reduces more deaths compared to the
contact-based strategy when the contagion rate is high and the vaccination
rate is low. However, as the vaccination rate becomes higher, the contact-
based strategy outperforms the fatality-based strategy.

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 using a metapopulation model, the density
of epidemic states in each subpopulation is tracked instead of tracking the
epidemic states of each individual. The density of states evolves due to the
interactions among subpopulations and interactions that occur within the same
subpopulation. Because metapopulation models have lower dimensions compared to
networks, they allow more exhaustive studies on the spread of epidemic
diseases. The epidemic equation for the SIRD model in the metapopulation model
is

$\displaystyle\frac{\partial}{\partial t}\rho_{\alpha}^{\mathrm{I}}(t)$
$\displaystyle=\eta\left(\rho_{\alpha}^{\mathrm{S}}(t)-v_{\alpha}\right)\sum_{\beta}M_{\alpha\beta}\rho_{\beta}^{\mathrm{I}}(t)-\mu\rho_{\alpha}^{\mathrm{I}}(t)\,,$
(4) $\displaystyle\frac{\partial}{\partial t}\rho_{\alpha}^{\mathrm{R}}(t)$
$\displaystyle=\left(1-\kappa_{\alpha}\right)\mu\rho_{\alpha}^{\mathrm{I}}(t)\,,$
(5) $\displaystyle\frac{\partial}{\partial t}\rho_{\alpha}^{\mathrm{D}}(t)$
$\displaystyle=\kappa_{\alpha}\mu\rho_{\alpha}^{\mathrm{I}}(t)\,,$ (6)

where $\rho_{\alpha}^{\mathrm{S}}$, $\rho_{\alpha}^{\mathrm{I}}$,
$\rho_{\alpha}^{\mathrm{R}}$, and $\rho_{\alpha}^{\mathrm{D}}$ are the
probabilities that an individual in group $\alpha$ is in the S, I, R, and D
state, respectively; $v_{\alpha}$ is the fraction of vaccinated individuals in
subpopulation $\alpha$; and $M_{\alpha\beta}$ is the contact matrix, which is
defined as the average contacts that an individual in group $\alpha$ has with
the individuals in group $\beta$.

Initially, an infinitesimal fraction, $n_{0}=10^{-8}$ of each group $\alpha$
of the population, is in the I state, and all the rest of the population,
$1-n_{0}$, is in the S state. As long as $n_{0}$ is small enough, the value of
$n_{0}$ and how these initially infected individuals are distributed among the
subpopulations do not affect the final states $\rho^{\mathrm{R}}$ and
$\rho^{\mathrm{D}}$. The differential equations are then solved by the fourth
order Runge-Kutta method [70, 71] until the total fraction of infected
individuals,
$\rho^{\mathrm{I}}=\sum_{\alpha}P_{\alpha}\rho_{\alpha}^{\mathrm{I}}$, becomes
less than a certain threshold, $10^{-12}$, and the epidemic process ends
($P_{\alpha}$ is the fraction of individuals in subpopulation $\alpha$.). We
then calculate the total fraction of the deceased population
$\rho^{\mathrm{D}}=\sum_{\alpha}P_{\alpha}\rho_{\alpha}^{\mathrm{D}}$.

We constructed a metapopulation model with heterogeneous contact and fatality
rates. The population has fatality rates $\kappa_{i}=$5%, 7.5%, 10%, 12.5%,
and 15% and relative contact rates $c_{j}=$0.5, 0.75, 1, 1.25, and 1.5. The
population is equally divided into 25 subpopulations according to the five
fatalities and five contact rates ($5\times 5=25$). The contact rate between
groups $(i,j)$ and $(i^{\prime},j^{\prime})$ is
$M_{iji^{\prime}j^{\prime}}=c_{j}c_{j}^{\prime}$.

We investigated the effectiveness of random, fatality-based, and contact-based
strategies for various levels of vaccine supply. In the random strategy, the
vaccine is randomly distributed and each subpopulation is uniformly
vaccinated. In the fatality-based strategy, the subpopulations are vaccinated
in descending order of fatality rates, and if two subpopulations have
identical fatality rates, the one with a higher contact rate is vaccinated. In
the contact-based method, an 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 $\alpha$ with unvaccinated individuals is

$\sum_{\beta}M_{\alpha\beta}\left(1-v_{\beta}\right)\,,$ (7)

and the value is recalculated at each iteration. This differs from the contact
rate of age group $\alpha$ with any individual in the population, which is
$\sum_{\beta}M_{\alpha\beta}$.

The mortality rate of the population when each vaccination strategy is
employed is illustrated in Figs. 1(a) and (b). When the contagion rate is low,
the contact-based strategy results in a lower mortality rate than the
fatality-based strategy regardless of the vaccination rate; however, for a
high contagion rate, there is a crossover between the strategies. The
fatality-based strategy more effectively reduces mortality when the
vaccination rate is low, but the contact-based strategy outperforms the
fatality-based strategy when the vaccination rate is high. If the vaccination
rate is sufficiently high, herd immunity is achieved regardless of the choice
of the vaccination strategy (fatality-based, contact-based, random, etc). The
difference between the mortality rates when fatality- and contact-based
strategies are employed is depicted in Fig. 1(c). The fatality-based strategy
is effective when the contagion rate is high and the vaccine supply is low.

Figure 2: (a, d) Fraction of recovered population, (b, e) average fatality of
the vaccinated population, and (c, f) average contact rate of the vaccinated
population of the synthetic metapopulation model with heterogeneous fatality
and contact rates. The recovery rate is normalized to $\mu=1$, and contagion
rates are (a–c) $\eta=0.05$ and (d–f) $\eta=0.4$. The increasing (decreasing)
curve, where the locally optimal vaccination strategies are found by
iteratively increasing (decreasing) the vaccination rate, is depicted as solid
red (dashed blue) lines. The globally optimal strategies are found by
simulated annealing (SA). There is no abrupt transition nor separation of the
curves for a low contagion rate, $\eta=0.05$, and the optimal vaccination
strategy prefers to vaccinate individuals with high contact rates regardless
of the vaccination rate. For a large contagion rate, $\eta=0.4$, there is an
abrupt transition in the globally optimal strategy. The separation of the
increasing and decreasing curve indicates the path-dependency of the
vaccination strategy which is analogous to hysteresis. For a small (0–0.56) or
large (0.67–1) vaccination rate, the increasing and decreasing curves
coincide; however, near the transition point, the increasing curve tends to
vaccinate individuals with high fatality rates (high-fatality strategy) and
the decreasing curve vaccinates individuals with high contact rates (high-
contact strategy). The High-fatality strategy results in a higher number of
recovered population than the high-contact strategy because the high-fatality
strategy aims to protect high-risk groups, while the high-contact strategy
aims to contain the infection.

Figure 3: The mortality rate of a mixed strategy that combines the high-
fatality and high-contact strategies. The contagion rate is $\eta=0.4$, and
vaccine supply is, from top to bottom, 58%, 60%, 62%, 64%, and 66%. There is a
barrier of mortality rate between the high-fatality and high-contact
strategies, which is the cause of the path-dependency and discontinuous
transition of the vaccination strategy.

The mortality rate of the high-contact strategy drops faster than the high-
fatality strategy as the level of vaccine supply increases, causing the
crossover between the two strategies.

###  III.2 Transition and path-dependency of the optimal vaccination strategy

In this section, we further investigate the vaccination rate dependency of the
vaccination strategy and demonstrate that the optimal vaccination strategy
undergoes a discontinuous transition. To find the globally optimal vaccination
strategy, we implement a modified version of the simulated annealing
technique. The simulated annealing is a probabilistic optimization algorithm
inspired by spin glass [72]. First, we start with a random vaccination
strategy with a given amount of vaccine supply. We set this strategy as a
provisional solution. We then 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{\left(-1/T\right)}$, where $T$ is the temperature of the
algorithm. In the beginning, the temperature is set at $T=2$. We iterate this
process $n_{\mathrm{iter}}=10^{6}$ times, while the temperature is dropped by
a factor of $f_{\mathrm{iter}}=1-2\times 10^{-5}$ at each step. The resulting
provisional solution is the optimal vaccination strategy, given that
$n_{\mathrm{iter}}$ is sufficiently large and $f_{\mathrm{iter}}$ is
sufficiently close to one. To find locally optimal solutions, we use the zero-
temperature simulated annealing, which is analogous to the gradient-descent
method. We perturb the provisional solution by decreasing the vaccination rate
of group $\alpha$ by $\delta v/P\left(\alpha\right)$ and increasing the
vaccination rate of group $\beta$ by $\delta v/P\left(\beta\right)$, where
$\delta v$ is a small number, and $P\left(\alpha\right)$ is the fraction of
the group $\alpha$ 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.

To investigate the path-dependency of the optimal vaccination strategy, we
iteratively increased and decreased the vaccination rate by a small amount,
while constantly calculating the locally optimal vaccination strategy in the
vicinity. To obtain the increasing curve, we first set the vaccination rate to
$\Delta v=0.01$ and find the optimal vaccination strategy
$v^{(I)}_{\beta}(\Delta v)$. We then increase the vaccination rate by $\Delta
v$ and find the locally optimal vaccination strategy $v^{(I)}_{\beta}(2\Delta
v)$ near the optimal strategy from the previous step. We repeat this process
until the vaccination rate reaches one. To obtain the decreasing curve, we
start from a vaccination rate of $1-\Delta v$ and repeat the process.

The results are illustrated in Fig. 2. Fraction of the recovered population
$\rho^{\mathrm{R}}$, average fatality, and contact rates of the vaccinated
individuals are depicted as characteristics of vaccination strategies. These
values are analogous to the order parameters of the phase transition in
thermal systems and the vaccination rate is the control parameter. The order
parameters of the two local mortality minima are depicted in the curves
similarly to the magnetization of the two free energy minima is depicted in
the hysteresis curve of the magnetic systems. The global mortality minimum
corresponds to the global free energy minimum where the system lies in the
Boltzmann distribution. For a low contagion rate, there is no abrupt
transition of the optimal vaccination strategy. For a high contagion rate, an
abrupt transition of the globally optimal vaccination strategy, which is
obtained by simulated annealing, discontinuously changes. Moreover, there is a
path-dependency in the vaccination strategy. When locally optimal vaccination
strategies are found by slowly increasing the vaccination rate from zero, the
strategies vaccinate individuals with high fatality rates in the middle region
(high-fatality strategy). If the strategies are found by slowly decreasing the
vaccination rate from one, they primarily vaccinate individuals with high
contact rates (high-contact strategy). A high-fatality strategy results in a
higher fraction of recovered individuals than the high-contact strategy, even
though the strategies’ mortality rates are similar or the same in the vicinity
of the transition point.

The path-dependency of this transition implies that a moderate strategy that
combines the high-fatality and high-contact strategy can be less effective
than either strategy. The vaccination rate of the moderate strategy is
$v_{\alpha}^{\mathrm{mod}}=rv_{\alpha}^{\mathrm{f}}+(1-r)v_{\alpha}^{\mathrm{c}}$,
where $v_{\alpha}^{\mathrm{f}}$ and $v_{\alpha}^{\mathrm{c}}$ are the
vaccination rates of subpopulation $\alpha$ in the high-fatality and high-
contact strategy, respectively. The performance of the moderate strategy for
various levels of vaccine supply is depicted in Fig. 3. There is a barrier of
mortality rate between the high-fatality and high-contact strategies, and the
moderate strategy is never more effective than both of the strategies, and in
some regions, it is less effective than either of the two strategies. Hence,
it is inadvisable to mix the two strategies or change from one to the other in
the middle. The mortality rate of the high-contact strategy ($r=0$) decreases
faster than the high-fatality strategy ($r=1$), which results in an abrupt
transition of the optimal vaccination strategy from a high-fatality to a high-
contact strategy.

###  III.3 Real-world epidemic diseases

Figure 4:  The age contact matrix of the (a) United Kingdom and (b) United
States, and (c) contact rate of each age group. The population is divided into
17 groups: aged 0–4, 5–9, …, 75–79, and above 80. The interaction strength
between groups of similar age is disproportionately higher, and the groups of
ages 10–29 show the highest number of contacts. (d) The fatality rates of TB
in the United Kingdom (UK) and COVID-19 in the United States (US). The
fatality rate is highly heterogeneous and monotonically increases with age,
with the sole exception of children below age five in TB. Therefore, the
senior population is primarily vaccinated by the fatality-based strategy, and
individuals of age 10–29 are primarily vaccinated by the contact-based
strategy.  Figure 5: Average age of vaccinated individuals for (a) TB, (b)
COVID-19, and (c) COVID-19 with reinfection and finite vaccine efficacy. The
contagion rates are $\eta=0.25$, and the recovery rate is normalized to
$\mu=1$. The total amount of vaccine is increased (solid red line) and
decreased (dashed blue line) while keeping the vaccination strategy at its
local optimum. For a small vaccination rate, the increasing curve and
decreasing curve coincide; however, for a sufficiently large vaccination rate,
the two curves show a largely unequal average age of the vaccinated
population. The increasing curve has a higher average age. The globally
optimal strategy, which is calculated by the simulated annealing, undergoes an
abrupt transition from the high-fatality to the high-contact strategies. The
separation of the increasing and decreasing curve indicates the path-
dependency of the vaccination strategy which is analogous to hysteresis.
Figure 6: The mortality rate of the mixed strategy of high-fatality and high-
contact strategies for (a) TB and (b) COVID-19. The contagion rate is
$\eta=0.25$, and the vaccine supply is, from top to bottom, 52%, 54%, 56%,
58%, and 60%. There is a barrier of mortality rate between the high-fatality
and high-contact strategy, which is the cause of the discontinuous transition
and path-dependency of the optimal vaccination strategy. The mortality rate of
the high-contact strategy drops faster than the high-fatality strategy as the
level of vaccine supply increases, and herd immunity is achieved at a lower
vaccine supply.

In this section, we show that the discontinuous transition and path-dependency
of the optimal vaccination strategy illustrated in the previous section occur
in the vaccination of actual epidemic diseases, such as tuberculosis (TB) and
COVID-19. Contact data between each age group in the various countries have
been studied utilizing surveys [73]. To model TB, we employed a contact matrix
calculated from the UK data along with the incidence risk ratio of TB in the
UK [74]. The population is divided into 17 groups: aged 0–4, 5–10, …, 75–79,
and above 80. For COVID-19, we use the contact matrix of the US and the age-
dependent IFR obtained from a meta-analysis of medical literature [50]. The
latter is calculated as

$\log_{10}{\rm IFR}=(-3.27\pm 0.07)+(0.0524\pm 0.0013)\,{\rm age}\,.$ (8)

The contact matrices and the fatality rates of the diseases are illustrated in
Fig. 4. Contact within a similar age group is disproportionately intense
(Figs. 4(a, b)), and contacts between teenagers exhibit the highest strength.
Fatality rates of the diseases monotonically increase with age, except for
children below age five for TB. As a result, the fatality-based strategy
primarily vaccinates the senior population, while the contact-based strategy
vaccinates the teenagers first.

The average age of the vaccinated individuals are presented for TB (Fig. 5(a))
and COVID-19 (Fig. 5(b)). Both figures exhibit the discontinuous transition
and path-dependency of the optimal vaccination strategy. The increasing curve
vaccinates the senior population more than the decreasing curve, which
corresponds to the phenomenon depicted in Fig. 2(e), where the increasing
curve has a greater preference to vaccinate individuals from groups with high
fatality rates than the decreasing curve. Also, there are barriers of
mortality rate between the high-fatality and high-contact strategies (Fig. 6),
suggesting that mixing the two strategies is ineffective. As the vaccination
rate increases, the mortality associated with the high-contact strategy
decreases faster than that of the high-fatality strategy to achieve herd
immunity at a lower vaccine supply.

###  III.4 Complex epidemic stages, vaccine breakthrough infection, and
reinfection

The previous results are obtained in simplified models. In this section, we
demonstrate that the more complicated behaviors of the actual diseases do not
significantly alter the findings of this research. First, in the real world,
actual infectious diseases progress in a series of epidemic stages, such as
the incubation period, prodromal period, and acute period. Each of these
stages has a distinct rate of spreading the disease. These stages have
complicated effects on the temporal dynamics of epidemics, but in this study,
only the fraction of population in each epidemic state (R and D) at the end of
the epidemic is relevant. In this sense, the complex stages of a disease can
be reduced to a simplified model. For instance, suppose there multiple
infectious stages $\mathrm{I}_{k}$ ($k=1,\cdots,K$) of a disease, each with
contagion rate $\eta_{k}$ and progression rate $\mu_{k}$ (i.e.,
$\mathrm{S}+\mathrm{I_{k}}\rightarrow\mathrm{I}_{1}+\mathrm{I}_{k}$ occurs
with rate $\eta_{k}$, $\mathrm{I}_{k}\rightarrow\mathrm{I}_{k+1}$ occurs with
rate $\mu_{k}$, and $\mathrm{I}_{K}\rightarrow\mathrm{R}$ occurs with rate
$\mu_{K}$). The total recovered population at the end of the epidemic disease
is then identical to that of the SIR model with contagion rate
$\eta^{*}=\sum_{k}\eta_{k}/\mu_{k}$ and $\mu=1$. To model an incubation stage,
we can set $\eta_{I}=0$ and $\mu_{I}=1/\tau_{I}$, where $\tau_{I}$ is the
incubation period of the disease.

Also, we assumed that vaccinated individuals never become infected even in
contact with infected individuals. However, people who are vaccinated still
can get infected by COVID-19. An infection of a vaccinated individual is
referred to as a vaccine breakthrough infection. To include vaccine
breakthrough infection, we can suppose a vaccine efficacy of $\theta<1$, and
the vaccinated individuals turn into the infected state at a rate of
$(1-\theta)\eta$ instead of $\eta$. Additionally, even when an infected
individual recovers and obtains immunity to the disease, there is a small
probability that the individual can be infected by the disease again. Such
reinfection can be modeled as some individuals losing immunity [75, 76, 77].
Hence, individuals in state R turn into S state at rate $\nu$. The typical
time for an individual to lose immunity is $1/\nu$. Individuals in the D state
remain in the D state.

We included vaccine breakthrough infections and reinfections to COVID-19 and
illustrated the average age of the vaccinated population in Fig. 5(c). The
vaccine efficacy is $\theta=0.9$, and the rate of immunity loss is $\nu=0.05$.
This means that typically the immunity of a vaccine is lost over a duration 20
times the average recovery time of the disease ($\sim 200$ days). The
discontinuous transition of the optimal strategy and path-dependency still
manifests themselves in the model with these modifications.

##  IV Conclusion

In summary, we employed the SIRD model to investigate the effectiveness of
vaccination strategies to minimize the mortality rate in a population with
heterogeneous fatality rates. We constructed a synthetic metapopulation model
with heterogeneous fatality and contact rates to investigate how the
effectiveness of vaccination strategies relates to the amount of vaccine
available. Vaccinating individuals with high fatality rates is effective when
the contagion rate is high and the vaccine supply is low. We found the
discontinuous transition and path-dependency, which is analogous to hysteresis
in statistical physics, of the optimal vaccination strategy. The path-
dependency of the vaccination strategy implies that combining high-fatality
and high-contact strategies is ineffective in reducing the mortality rate of
the epidemic disease. We also demonstrated that such phenomena occur in real-
world epidemic diseases, such as TB and COVID-19. These conclusions are valid
even when complex stages of a disease, vaccine breakthrough infection, and
reinfection are considered.

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 were treated equally. However, the contagion rate of disease
among 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 are collected and implemented, the relative
strength of such interactions can 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 discontinuous
transition and path-dependency of the optimal vaccination strategies occur in
various epidemic models with a wide range of parameters, the conclusions of
this research should still be valid.

###### Acknowledgements.

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

## Data Availability Statement

The data that support the findings of this study are openly available in
https://github.com/mobs-lab/mixing-patterns [73] (contact matrices), Ref. [74]
(age-specific fatality of TB), and Ref. [50] (age-specific fatality of
COVID-19).

## References

  * Dorogovtsev  _et al._ [2008] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008).
  * Zimmermann  _et al._ [2004] M. G. Zimmermann, V. M. Eguíluz, and M. San Miguel, Phys. Rev. E 69, 065102 (2004).
  * Dorogovtsev and Mendes [2002] S. N. Dorogovtsev and J. F. Mendes, Adv. Phys. 51, 1079 (2002).
  * Klemm  _et al._ [2003] K. Klemm, V. M. Eguíluz, R. Toral, and M. San Miguel, Phys. Rev. E 67, 026120 (2003).
  * Nekovee  _et al._ [2007] M. Nekovee, Y. Moreno, G. Bianconi, and M. Marsili, Physica A 374, 457 (2007).
  * Watts  _et al._ [2005] D. J. Watts, R. Muhamad, D. C. Medina, and P. S. Dodds, Proc. Natl. Acad. Sci. U.S.A. 102, 11157 (2005).
  * Colizza and Vespignani [2007] V. Colizza and A. Vespignani, Phys. Rev. Lett. 99, 1 (2007).
  * Colizza and Vespignani [2008] V. Colizza and A. Vespignani, J. Theor. Biol. 251, 450 (2008).
  * Lloyd and Jansen [2004] A. L. Lloyd and V. A. Jansen, Math. Biosci. 188, 1 (2004).
  * Masuda [2010] N. Masuda, New J. Phys. 12, 10.1088/1367-2630/12/9/093009 (2010).
  * Hoang and Antoncic [2003] H. Hoang and B. Antoncic, J. Bus. Ventur. 18, 165 (2003).
  * Newman [2003] M. E. Newman, SIAM Rev. 45, 167 (2003).
  * Boccaletti  _et al._ [2006] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Phys. Rep. 424, 175 (2006).
  * Acemoǧlu  _et al._ [2013] D. Acemoǧlu, G. Como, F. Fagnani, and A. Ozdaglar, Math. Oper. Res. 38, 1 (2013).
  * Grabowski and Kosiński [2006] A. Grabowski and R. A. Kosiński, Physica A 361, 651 (2006).
  * Watts and Dodds [2007] D. J. Watts and P. S. Dodds, J. Consum. Res. 34, 441 (2007).
  * Moreno  _et al._ [2002] Y. Moreno, R. Pastor-Satorras, and A. Vespignani, Eur. Phys. J. B 26, 521 (2002).
  * Ji and Jiang [2014] C. Ji and D. Jiang, Appl. Math. Model. 38, 5067 (2014).
  * Pastor-Satorras and Vespignani [2001a] R. Pastor-Satorras and A. Vespignani, Phys. Rev. E 63, 66117 (2001a).
  * Mata and Ferreira [2013] A. S. Mata and S. C. Ferreira, EPL 103, 48003 (2013).
  * de Arruda  _et al._ [2020] G. F. de Arruda, G. Petri, F. A. Rodrigues, and Y. Moreno, Phys. Rev. Res. 2, 13046 (2020).
  * Katona  _et al._ [2011] Z. Katona, P. P. Zubcsek, and M. Sarvary, J. Mark. Res. 48, 425 (2011).
  * Rogers [2004] E. M. Rogers, J. Health Commun. 9, 13 (2004).
  * Jhun  _et al._ [2019] B. Jhun, M. Jo, and B. Kahng, J. Stat. Mech. 2019, 123207 (2019).
  * Jhun [2021] B. Jhun, Phys. Rev. Res. 3, 033282 (2021), arXiv:2108.04568 .
  * Jia and Benson [2019] J. Jia and A. R. Benson, in  _Proc. Twelfth ACM Int. Conf. Web Search Data Min._ (ACM, New York, NY, USA, 2019) pp. 366–374.
  * Leskovec  _et al._ [2007] J. Leskovec, J. Kleinberg, and C. Faloutsos, ACM Trans. Knowl. Discov. Data 1, 10.1145/1217299.1217301 (2007).
  * Fowler [2006] J. H. Fowler, Polit. Anal. 14, 456 (2006).
  * Pastor-Satorras and Vespignani [2001b] R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001b).
  * Choi  _et al._ [2017] W. Choi, D. Lee, and B. Kahng, Phys. Rev. E 95, 22304 (2017).
  * Lee  _et al._ [2017] D. Lee, W. Choi, J. Kertész, and B. Kahng, Sci. Rep. 7, 1 (2017).
  * Lee  _et al._ [2016] D. Lee, S. Choi, M. Stippinger, J. Kertész, and B. Kahng, Phys. Rev. E 93, 42109 (2016).
  * Pastor-Satorras and Vespignani [2002] R. Pastor-Satorras and A. Vespignani, Phys. Rev. E 65, 1 (2002).
  * Hébert-Dufresne  _et al._ [2013] L. Hébert-Dufresne, A. Allard, J.-G. Young, and L. J. Dubé, Sci. Rep. 3, 2171 (2013).
  * Cohen  _et al._ [2003] R. Cohen, S. Havlin, and D. Ben-Avraham, Phys. Rev. Lett. 91, 2 (2003).
  * Ghalmane  _et al._ [2019] Z. Ghalmane, M. E. Hassouni, and H. Cherifi, Soc. Netw. Anal. Min. 9, 1 (2019).
  * Osat  _et al._ [2017] S. Osat, A. Faqeeh, and F. Radicchi, Nat. Commun. 8, 1 (2017).
  * Schneider  _et al._ [2011] C. M. Schneider, T. Mihaljev, S. Havlin, and H. J. Herrmann, Phys. Rev. E 84, 061911 (2011).
  * Costa and Ferreira [2020] G. S. Costa and S. C. Ferreira, Phys. Rev. E 101, 022311 (2020).
  * Yan  _et al._ [2015] S. Yan, S. Tang, W. Fang, S. Pei, and Z. Zheng, J. Stat. Mech. 2015, 10.1088/1742-5468/2015/08/P08010 (2015).
  * Clusella  _et al._ [2016] P. Clusella, P. Grassberger, F. J. Pérez-Reche, and A. Politi, Phys. Rev. Lett. 117, 1 (2016).
  * Masuda [2009] N. Masuda, New J. Phys. 11, 0 (2009).
  * Van Mieghem  _et al._ [2011] P. Van Mieghem, D. Stevanović, F. Kuipers, C. Li, R. Van De Bovenkamp, D. Liu, and H. Wang, Phys. Rev. E 84, 1 (2011).
  * Matamalas  _et al._ [2018] J. T. Matamalas, A. Arenas, and S. Gómez, Sci. Adv. 4, eaau4212 (2018).
  * Dong  _et al._ [2020] Z. Dong, Y. Chen, T. S. Tricco, C. Li, and T. Hu, in  _2020 IEEE Intl Conf Parallel Distrib. Process. with Appl. Big Data Cloud Comput. Sustain. Comput. Commun. Soc. Comput. Netw._ (IEEE, 2020) pp. 845–851.
  * Chen  _et al._ [2008] Y. Chen, G. Paul, S. Havlin, F. Liljeros, and H. E. Stanley, Phys. Rev. Lett. 101, 2 (2008).
  * Madar  _et al._ [2004] N. Madar, T. Kalisky, R. Cohen, D. Ben-Avraham, and S. Havlin, Eur. Phys. J. B 38, 269 (2004).
  * Tatar  _et al._ [2021] M. Tatar, J. M. Shoorekchali, M. R. Faraji, and F. A. Wilson, J. Glob. Health 11, 03086 (2021).
  * Li  _et al._ [2020] W. Li, R. Thomas, H. El-Askary, T. Piechota, D. Struppa, and K. A. Abdel Ghaffar, Earth Syst. Environ. 4, 513 (2020).
  * Levin  _et al._ [2020] A. T. Levin, W. P. Hanage, N. Owusu-Boaitey, K. B. Cochran, S. P. Walsh, and G. Meyerowitz-Katz, Eur. J. Epidemiol. 35, 1123 (2020).
  * Manuel  _et al._ [2020] B. Manuel, K. Richard, T.-s. Sarah, H. H. H, W. A. F, and N. R. A, Swiss Med Wkly 150, 2019 (2020).
  * Barone-Adesi  _et al._ [2020] F. Barone-Adesi, L. Ragazzoni, and M. Schmid, Disaster Med. Public Health Prep. 14, E1 (2020).
  * Kim  _et al._ [2020] D. H. Kim, Y. J. Choe, and J. Y. Jeong, J. Korean Med. Sci. 35, 1 (2020).
  * Shim [2021] E. Shim, Int. J. Environ. Res. Public Health 18, 10.3390/ijerph18105053 (2021).
  * Bhatt  _et al._ [2021] T. Bhatt, V. Kumar, S. Pande, R. Malik, A. Khamparia, and D. Gupta, Stud. Comput. Intell. 924, 25 (2021).
  * Scarpino and Petri [2019] S. V. Scarpino and G. Petri, Nat. Commun. 10, 10.1038/s41467-019-08616-0 (2019).
  * Cai  _et al._ [2015] W. Cai, L. Chen, F. Ghanbarnejad, and P. Grassberger, Nat. Phys. 11, 936 (2015).
  * Chen  _et al._ [2013] L. Chen, F. Ghanbarnejad, W. Cai, and P. Grassberger, EPL 104, 10.1209/0295-5075/104/50001 (2013).
  * Cho  _et al._ [2009] Y. S. Cho, J. S. Kim, J. Park, B. Kahng, and D. Kim, Phys. Rev. Lett. 103, 135702 (2009).
  * Bianconi and Krapivsky [2020] G. Bianconi and P. L. Krapivsky, Phys. Rev. E 102, 1 (2020).
  * Choi  _et al._ [2022] K. Choi, H. Choi, and B. Kahng, Chaos, Solitons & Fractals 157, 111904 (2022), arXiv:2010.07157 .
  * Scarselli  _et al._ [2021] D. Scarselli, N. B. Budanur, M. Timme, and B. Hof, Nat. Commun. 12, 10.1038/s41467-021-22725-9 (2021).
  * Li  _et al._ [2022] W. Li, X. Xue, L. Pan, T. Lin, and W. Wang, Appl. Math. Comput. 412, 126595 (2022), arXiv:2105.08234 .
  * St-Onge  _et al._ [2022] G. St-Onge, I. Iacopini, V. Latora, A. Barrat, G. Petri, A. Allard, and L. Hébert-Dufresne, Commun. Phys. 5, 25 (2022).
  * Liu  _et al._ [2022] Y. Liu, B. Liu, Y. Deng, and J. Liu, Complexity 2022, 1 (2022).
  * Wang and Jia [2019] P. Wang and J. Jia, Adv. Differ. Equations 2019, 433 (2019).
  * Sen and Sen [2021] D. Sen and D. Sen, Ind. Eng. Chem. Res. 60, 4251 (2021).
  * Yuan  _et al._ [2020] M. Yuan, W. Yin, Z. Tao, W. Tan, and Y. Hu, PLoS One 15, 1 (2020).
  * Anastassopoulou  _et al._ [2020] C. Anastassopoulou, L. Russo, A. Tsakris, and C. Siettos, PLoS One 15, 1 (2020).
  * Runge [1895] C. Runge, Math. Ann. 46, 167 (1895).
  * Kutta [1901] W. Kutta, Zeit. Math. Phys. 46, 435 (1901).
  * Kirkpatrick  _et al._ [1987] S. Kirkpatrick, C. Gelatt, and M. Vecchi, in  _Readings Comput. Vis._, 4598 (Elsevier, Berlin, Heidelberg, 1987) pp. 606–615.
  * Mistry  _et al._ [2021] D. Mistry, M. Litvinova, A. Pastore y Piontti, M. Chinazzi, L. Fumanelli, M. F. Gomes, S. A. Haque, Q. H. Liu, K. Mu, X. Xiong, M. E. Halloran, I. M. Longini, S. Merler, M. Ajelli, and A. Vespignani, Nat. Commun. 12, 1 (2021).
  * Pedrazzoli  _et al._ [2019] D. Pedrazzoli, K. Kranzer, H. L. Thomas, and M. K. Lalor, ERJ Open Research 5, 10.1183/23120541.00125-2019 (2019).
  * Okuwa  _et al._ [2019] K. Okuwa, H. Inaba, and T. Kuniya, Math. Biosci. Eng 16, 6071 (2019).
  * Saif [2019] M. A. Saif, Physica A: Statistical Mechanics and its Applications 535, 122251 (2019).
  * Salman  _et al._ [2021] A. M. Salman, I. Ahmed, M. H. Mohd, M. S. Jamiluddin, and M. A. Dheyab, Computers in Biology and Medicine 133, 104372 (2021).