Pair approximation model for the vaccination game: predicting the dynamic process of epidemic spread and individual actions against contagion

We successfully establish a theoretical framework of pairwise approximation for the vaccination game in which both the dynamic process of epidemic spread and individual actions in helping prevent social behaviours are quantitatively evaluated. In contrast with mean-field approximation, our model captures higher-order effects from neighbours by using an underlying network that shows how the disease spreads and how individual decisions evolve over time. This model considers not only imperfect vaccination but also intermediate protective measures other than vaccines. Our analytical predictions are validated by multi-agent simulation results that estimate random regular graphs at varying degrees.


Introduction
In today's globalized world, all kinds of infectious diseases, such as measles, smallpox, pertussis, flu and Ebola, can spread rapidly around the world and negatively affect the lives of people. Vaccination is one of the most effective measures for preventing the transmission of infectious diseases and for reducing morbidity and mortality rates [1,2]. However, voluntary vaccination policies create a social dilemma [3,4], namely, the 'vaccination dilemma' or 'paradox of epidemiology'. As a result, individuals often hesitate to get themselves vaccinated. When herd immunity is achieved in a community with increased vaccination coverage (VC), non-vaccinators can avoid infections and vaccination costs. However, herd immunity is inevitably destroyed because nonvaccinators completely depend on the efforts of vaccinators. Therefore, the actual VC will not be able to satisfy demand.
To model this vaccination dilemma, researchers have studied vaccination games, which can predict the dynamics of (i) epidemic spread in a complex social network and (ii) decision making on whether to undergo vaccination depending on the status of the epidemic. Epidemic dynamics are predicted using mathematical epidemiological models, such as the susceptibleinfected-recovered (SIR) model, whereas decision making is modelled on evolutionary game theory. Bauch et al. [5] and Fu et al. [6] published pioneering works on the vaccination game. To quantitatively investigate the multiple effects on vaccination behaviour, a significant number of researchers have studied various frameworks [7][8][9][10][11][12][13]. Most studies on the vaccination game have relied on multi-agent simulation (MAS), which allows for a more flexible and realistic modelling approach. In addition to the MAS approach, Fu et al. [6] proposed a mathematical framework for a mixed-population vaccination game that assumes perfect vaccination. A theoretical approach based on a set of ordinary differential equations (ODEs) can be a powerful tool for explicitly demonstrating the dynamics of both epidemic spread and human decision making. In most recent study, Steinegger et al. [14] provided a theoretical framework that allows to incorporate the structure of the network in this kind of disease-driven dilemmas using Markovian equations.
Most studies on the vaccination game have assumed that vaccinations provide perfect immunity to each vaccinator. In reality, vaccinations can only impart partial protection against many infectious diseases, such as measles, influenza, malaria and HIV. In addition to vaccinations, there are other protective measures from infectious diseases, such as mask-wearing, gargling and hand washing, which are called intermediate protective measures. Although these protective measures come at a more reasonable cost than vaccinations, they cannot block the transmission of infections to the body as effectively as vaccinations. Therefore, the stochastic effects of imperfect vaccination and intermediate protective measures need to be considered. On the basis of this background, Cardillo et al. [15] analysed the effects of imperfect vaccination on immunization behaviour in Erdős-Rényi random graph (ER-RG) and Barabasi-Albert scale-free (BA-SF) networks by using the MAS approach. Iwamura et al. [16] and Ida and Tanimoto [17] used MAS to investigate the effect of intermediate protective measures on square lattice and BA-SF networks. Wu et al. [18] and Kuga and Tanimoto [19] developed a new mathematical framework of the vaccination game that considered imperfect vaccination in an infinite and well-mixed population corresponding with a perfect graph by using mean-field approximation. Furthermore, they proposed a mathematical model for reproducing the vaccination game in complex networks corresponding to ER-RG and BA-SF networks by using degree distribution [20]. Several studies have recently used the same concept and investigated the multiple effects of imperfect vaccination and other parameters on vaccination behaviour [21][22][23][24][25][26][27][28][29][30][31][32]. Therefore, the modelling and analysis of the vaccination game can be enriched by the abundant quantity and quality of studies that have followed a MAS or theoretical approach.
While highlighting the theoretical approach, it should be noted that most mathematical studies on the vaccination game have focused on the so-called mean-field approximation, which is a single-order approximation for predicting the dynamics of epidemics and decision making, despite being extended to some heterogeneous topologies, such as infinite-sized graphs with degree distributions obeying the power law or Poisson distribution [20,30]. However, mean-field approximation has substantial limitations because it ignores higher order effects and assumes that one's neighbour can always be defined as a socially averaged hypothetical neighbour, which can never be true in a real-world, complex social network. This defect means that predictions by mean-field approximation are slightly different from those by MAS, which can capture predictions to a reasonable extent. As an alternative to mean-field approximation, the current study uses pair approximation, which has been widely applied in other fields. For example, in the spatial version of the Prisoner's Dilemma, which is an archetype model for quantifying network

Model set-up and results (a) Epidemic models
The pair approximation SIR model [35,37] was used to accurately model epidemic spread. We integrated the pair approximation SIR model with the concept of effectiveness and efficiency models, which can reproduce imperfect vaccination and intermediate protective measures [19,20]. The electronic supplementary material describes the fundamental formulation when the pair approximation is applied to SIR dynamics.

(i) Effectiveness of vaccination
In the effectiveness model, a vaccinated population is separated into two classes: individuals with perfect immunity and non-immune individuals. The population is subdivided into the following: non-vaccinated susceptible individual S N , non-vaccinated infected individual I N , non-vaccinated recovered individual R N , vaccinated susceptible individual S V , vaccinated infected individual I V , vaccinated recovered individual R V and vaccinated individual with perfect immunity P V . Let the effectiveness of the vaccination and the vaccination coverage be e (0 ≤ e ≤ 1) and x, respectively. The fraction of vaccinated individuals with perfect immunity [P V ](t) must be ex, On the basis of the above assumptions, the dynamics of the SIR/V model with imperfect vaccination can be described by the following ODEs: The above set of dynamic equations should be assumed to have the following set of initial conditions: Here, α is the vaccinator-non-vaccinator connection coefficient, which can be said dissortativity observed at the initial moment of every season (time-evolved in repeating seasons), and is less than x and 1 − x. If the homogeneous distribution of vaccinator and non-vaccinated was assumed, α = x(1 − x).
The following constraints are required: Solving the above set of equations, we obtain the following equations (the detailed mathematical manipulations are shown in the electronic supplementary material): , (2.20) At the steady state (t → ∞), the constraints in equation (2.14) to (2.16) can be rewritten as follows: and (2.28) 1/Q and taking into account the definition of μ = (Q − 1)/Q, we can write equations (2.27) and (2.28) as follows: The final fractions are expressed as follows: The other final pairs are hypothetically calculated as follows: The above set of equations are assumed to have the following initial conditions: The following constraints are required: Solving the above set of equations, we obtain the following equations (the detailed mathematical manipulations are shown in the electronic supplementary material): At the steady state (t → ∞), the constraints in equations (2.54) to (2.56) can be rewritten as follows: The final fractions are expressed as follows: The other final pairs are hypothetically calculated as follows:  Table 1 summarizes the payoff on the basis of whether the approach was V (either vaccination or intermediate protective measure) or NV and whether the individual is healthy or infected. On the basis of the expected payoff, the average social payoff π is formulated as follows: (Effectiveness model) (2.80) (Efficiency model) payoff in the previous epidemic season. In this study, we consider the strategy-updating approach proposed by Fu et al. [6], which is called 'individual-based risk assessment' (IB-RA).
In IB-RA, for two-strategy and two-player (2 × 2) games, an individual can update their strategy by randomly selecting neighbouring individuals stochastically and copying them. The individual decides whether to adopt that neighbour's strategy by using pairwise Fermi updating [6], i.e. individual i adopts the selected neighbour j's strategy with a probability: where s i is the strategy of i, π i is i's payoff in the previous season and parameter κ > 0 characterizes the strength of selection (the sensitivity of individuals to differences in their payoffs); a smaller κ indicates that an individual is more sensitive to a payoff difference. We set κ = 0.1 on the basis of [6]. The transition probability affecting the dynamics of x, which should be considered in light of the IB-RA rule, is covered by one of the following eight cases: (2.83) After the end of every epidemic season, the vaccination coverage x will increase or decrease.
To quantify this evolutionary process, we obtain the following equation for the dynamic system: As a result of changing the strategy, the vaccinator-non-vaccinator connection coefficient α will also change. To quantify this evolutionary process, we obtain the following equation for the dynamic system: Effectiveness model Here, P(A → VorNV|AB) is a transition probability that the focal A of pair AB change to the opposite strategy (vaccinator or non-vaccinator). All transition probabilities are described in the electronic supplementary material.  All of the above dynamic equations can be solved numerically. Therefore, the final result is affected by a two-stage process: single-season SIR dynamics and the strategy adaptation process. We rely on an explicit scheme, and we evaluate how this specific dynamic system evolves.   Therefore, we can observe the final epidemic size, vaccination coverage and average social payoff in social equilibrium.  an increase in degree, the present model's prediction gets close to the prediction of the meanfield approximation. A random regular graph with the infinite degree is quite possible because it literally means 'well mixed.' Keeping this in mind, we were able to see how meaningful the present model is over conventional analytical approaches that rely on mean-field approximation. In all figures, the red colour in final epidemic size, blue colour in vaccination coverage and red colour in average social payoff indicate a pandemic state, in which most individuals rely on neither vaccination nor intermediate protective measures. Therefore, an almost full-scale spread of infection is inevitable. Generally speaking, these regions emerge when the reliability of vaccinations or intermediate protective measure is low or when the cost is high. The border between each of these monotone regions and the remaining area represents a threshold that suggests a combination of critical effectiveness (efficiency) and critical cost to appropriately control the epidemic spread. Therefore, under a voluntary vaccination policy, we can confirm how the relationship between the reliability of protective measures and the cost of protection is an important factor in controlling the damage caused by the epidemic spread. Figure 3 shows the relationship between cost and critical effectiveness (efficiency) under which no individual takes vaccination (or intermediate protective measures). It clearly suggests that a higher cost requires a higher effectiveness (efficiency) to enable individuals to commit because individuals have no incentive to take protective measures unless there is a sufficiently high counterbalancing effect to the high cost. An observation of the differences in the degree Q shows that a lower degree requires greater effectiveness (efficiency), i.e. lower level social networks are more robust against  disease spreading than those at a higher level. This ironically results in a reduced incentive for individuals to vaccinate (or take intermediate protective measures) because of the temptation to take a free ride on herd immunity. As shown in figure 1 and 2, a lower effectiveness (or efficiency) corresponds with a higher vaccination coverage as long as the costs imposed are acceptable. This ironic situation can be explained as follows: when a protective measure is less reliable, it will create more uncertainty and fear, and more individuals will commit to the measure. However, even if a large proportion of individuals take protective measures, the epidemic cannot be eradicated if the vaccination is unreliable.

Discussion
A comparison of the effectiveness and efficiency models in figures 1 and 2 show that the latter has a wider pandemic phase at first glance. This implies that intermediate protective measures with certain η are less effective for suppressing an epidemic than imperfect vaccination, and this finding is consistent with our previous results based on mean-field approximation. Figure 4 shows the isograms of the relative vaccination cost C r drawn on a two-dimensional plane of vaccination coverage and effectiveness (efficiency). The grey arrows indicate the direction of the rising slope of C r , thus indicating that vaccination coverage (the rate of intermediate protective measures) gradually increases with decreasing effectiveness (efficiency) from e = 1.0 (η = 1.0) and dramatically drop down to zero at the critical effectiveness (efficiency), when the relative vaccination costs are low. This finding confirms what we have discussed above and is consistent with the reports by Wu et al. [18] and Chen et al. [21]. On the other hand, in case of high  To validate the theoretical results, a series of numerical simulations based on the MAS approach was performed. Following the procedure in previous studies [8], we set β as the point at which the final epidemic size exceeds the predefined threshold (i.e. 0.9) without any vaccinated individuals. Accordingly, we set β = 0.72, 0.38 and 0.14 for random regular graph with degree Q = 3, 4 and 8, respectively. We set γ = 1/3. A typical flu is assumed to determine these disease parameters. The results shown below were obtained by a collective average of 100 independent realizations starting from various initial conditions. The results for the effectiveness and efficiency models are shown in figures 5 and 6. Generally, all results are consistent with figures 1 and 2, although there are subtle discrepancies arising from the fact that the above simulation assumed a finite population size of N = 10 4 .
To characterize the effect of infection rates on vaccination behaviour, final epidemic size and vaccination coverage were shown as a function of the inverse of the effective recovery rate 1/r and the protection quality: (i) effectiveness e and (ii) efficiency η in figure 7. We fixed the relative vaccination cost C r = 0.   (1 − x)). We assumed perfect vaccination. (Online version in colour.) individuals do not vaccinate until the relative recovery rate r is less than the critical relative recovery rate r c , as expressed in Eq. (A-19) of the electronic supplemental material, because the epidemic does not spread. Furthermore, vaccination coverage increases with an increase in Q, both in terms of effectiveness and efficiency. This implies that a high number of degrees promote not only the epidemic spread but also the (cooperative) vaccination behaviour. The critical value of effectiveness, below which no individual will be vaccinated, does not show sensitivity to the relative cost, which is expressed as e = 0.35 in this case. On the other hand, the critical value of efficiency, below which no individual will be vaccinated, increases monotonically with an increase in the relative recovery rate. This result is consistent with results by Cardillo et al. [14].
In general, most of theoretical vaccination game model using mean-field approximation ignored the local correlations between vaccinators and non-vaccinators and assumed a homogeneous distribution of vaccinator and non-vaccinator. On the other hand, pair approximation can capture the evolution of heterogeneous distribution of vaccinator and non-vaccinator caused by strategy adaptation. In the model formulation, we introduced new parameter; α, meaning disassortativity at the beginning of every season, that can be time-evolved according to equations (2.86) and (2.87) when we digress from the mean-field approximation mentioned above. Let us name this setting 'heterogeneous' case. Contrariwise, when we follow to the conventional idea; taking mean-field approximation, α is always frozen at α = x(1 -x), which is called 'homogeneous' case. As shown in figure 8, when the heterogeneous setting was imposed, vaccination coverage is relatively high and final epidemic size is relatively low compared to the homogeneous one (α = x(1 − x)). Regardless of C r , we confirmed that α is always lower than that the homogeneous case; x(1 -x), because of the existence of the cluster of vaccinators and nonvaccinators. The existence of non-vaccinators cluster with high-frequency pushes up the risk of the disease spreading shared amid non-vaccinators, which makes defectors commit vaccination. This subsequently enables less final epidemic size when the heterogeneous case is presumed.

Conclusions
We presented a more precise theoretical framework of the vaccination game that considers imperfect vaccination and intermediate protective measures. The greatest advantage of our model is that it relies on a pair approximation epidemic model instead of mean-field approximation. The exact mathematical formulae for both dynamic processes, namely, epidemic spreading and strategy updating, are explicitly discussed. When solving the ODEs for the constructed pair approximation epidemic model, the critical vaccination coverage and the final fractions for each individual were derived. The results show that it is important to consider the degree effect in pairwise or mean-field approximation, particularly when a lesser degree is presumed. In addition, the effect of imperfect vaccination and intermediate protective measure strongly affects vaccination behaviour and final epidemic size. We validated our framework by comparing it with the MAS approach.
Data accessibility. All data have been generated from a series of numerical simulations. We have provided the source code (c++) for the numerical simulation in the electronic supplementary material.