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

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

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

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.

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.

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

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.

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.