Voluntary vaccination may not stop monkeypox A game-theoretic model

Authors: Ian B. Augsburger, Grace K. Galanthay, Jacob H. Tarosky, Jan Rychtář, Dewey Taylor

Abstract

Monkeypox (MPX) is a viral zoonotic disease that was endemic to Central and West Africa. However, during the first half of 2022, MPX spread to almost 60 countries all over the world. Smallpox vaccines are about 85% effective in preventing MPX infections. Our objective is to determine whether the vaccines should be mandated or whether voluntary use of the vaccine could be enough to stop the MPX outbreak. We incorporate a standard SVEIR compartmental model of MPX transmission into a game-theoretical framework. We study a vaccination game in which individuals decide whether or not to vaccinate by assessing their benefits and costs. We solve the game for Nash equilibria, i.e., the vaccination rates the individuals would likely adopt without any outside intervention. We show that, without vaccination, MPX can become endemic in previously non-endemic regions, including the United States. We also show that to “not vaccinate” is often an optimal solution from the individual’s perspective. Moreover, we demonstrate that, for some parameter values, there are multiple equilibria of the vaccination game, and they exhibit a backward bifurcation. Thus, without centrally mandated minimal vaccination rates, the population could easily revert to no vaccination scenario.

1 Introduction

1.1 Monkeypox

Monkeypox (MPX) is a viral zoonotic disease endemic to Central and West Africa [1]. The MPX cases suffer from mild symptoms such as headaches, fevers, rashes, lesions in their mouth and on their body [2], although there may be other and potentially severe complications such as blindness [3] or death [4]; see [5] for a comprehensive review.

It has recently garnered much public attention due to its 2022 outbreaks. From January 1 to July 4, 2022, 6027 laboratory confirmed MPX cases were reported to WHO from 59 countries, most of which are considered non-endemic to MPX [6]. Men who have sex with men (MSM) with new or multiple partners are amongst the most affected [7].

MPX is caused by a virus similar to smallpox virus and the smallpox vaccines provide about 85% protection from MPX [8]. For the current outbreak, [9] recommends (a) post-exposure prophylaxis (PEP) with an appropriate second- or third-generation vaccine for contacts of cases, and (b) pre-exposure prophylaxis (PrEP) for people at risk. However, as the MPX epidemic continues to unfold, there are calls to use the smallpox vaccine as PrEP in MSM at high risk of monkeypox virus exposure as it may also reduce transmission into the general population [7].

1.2 Mathematical models of MPX

Mathematical modeling is now a standard tool for disease prevention and elimination efforts [10, 11]. There are very few mathematical models specific for MPX; yet, in recent years, even before the 2022 outbreak, the modelling activity has been picking up. In [12], the authors developed the first model to represent MPX and other pox-like infections. The model was later extended to include the coexistence of HIV and MPX [13]. Culling as a means to prevent MPX was investigated in [14]. In [15] and [16], the authors studied transmission dynamics with treatment and vaccination; in [17], the authors developed a model for diagnosing MPX and in [18], the authors performed stability analysis for equilibria of their ODE system. Other models include [19] where the authors studied quarantining and public education, [20] and [21] which concerns human-to-human transmissions. The global and local asymptotic stability and transmission dynamics were explored in [22]. The impact of smallpox vaccines on MPX epidemics was investigated in [23]. The fractional order ODE models were developed in [24] and [25]. A model incorporating sexual behavior dynamics and stratifying the population into high- and low-risk groups was developed in [26]. Stochastic models and individual based simulations for the current outbreak are also being developed [27–30].

The potential for the disease spread in the population can be measured by the basic reproduction number, R0, the number of secondary infections from a single infected individual in a susceptible population [31]. If interventions in disease transmission are implemented, the number of secondary infections from a single infected individual in an otherwise healthy population is called the effective reproduction number, Re [32]. When R0, or more generally Re is less than 1, then the disease cannot spread in the population [33].

1.3 Game theory

In the context of this paper, a game is a mathematical model of a situation where several individuals interact (directly or indirectly) with one another and where each individual acts in its own interest [34]. The game theory has a rich history; the modern treatments build on the ideas of John von Neumann [35] and John Nash [36]; see for example [37–40] or, for more biologically oriented applications, [41–45].

One way to define the solution of a non-cooperative game is the so-called Nash equilibrium (NE). In the NE, each player is assumed to know the equilibrium strategies of the other players, and no one has any incentive to change only their own strategy [36]. In this paper, we will focus on this concept of NE, but we note that there are other possible approaches, involving for example bounded rationality [46, 47] or quantal response equilibria [48, 49].

Game theoretical models proved themselves useful in studying complex epidemiological scenarios in which self-interested individuals take actions based on the decisions of the rest of the population [50, 51]; and, as argued in [52], by incorporating human behavior into the epidemiological models one can get better insight and predictions.

1.4 Vaccination games

Vaccination games are games in which individuals decide whether to vaccinate or not. They are a class of public goods games [53] because vaccination produces public goods (herd immunity against a disease) that have the following two main non-rivalry, i.e., consumption of a good by one person does not affect other individuals, and non-exclusion of consumption, i.e., it is impossible to restrict the benefits to certain individuals [54]. Vaccination is prone to free-riding; the “free-riders” avoid the costs associated with vaccination while benefiting from vaccines taken by others [55]. People balance the perceived costs against the vaccine’s effectiveness [56] and it is well known that individuals act in a way that maximizes their self-interests, rather than the interests of the entire group [57].

The vaccination games have been applied to protection strategies to control diseases such as smallpox [58, 59], chickenpox [60], polio [61], influenza [62], Ebola [63], COVID-19 [64–66], chikungunya [67], Hepatitis B [68], lymphatic filiarisis [69]. [70] used the game theoretic framework to assess vaccination strategies with the presence of animal reservoirs of infection.

1.5 Content of the current paper

In the current paper, we extend the analysis from [70] by explicitly considering the MPX vaccine to be imperfect and allowing for the possibility of infections after vaccination. We focus on human-to-human transmission as done in [21]. We study Nash equilibria, the solutions of the vaccination game in which susceptible individuals decide whether or not to vaccinate against MPX. We calibrate our model based on historical data about MPX and we also use data from the 2022 outbreaks. The analysis reveals a possibility of multiple Nash equilibria and the existence of backward bifurcation. We perform sensitivity analysis and also study a hypothetical scenario under which the MPX transmission rate is higher than generally assumed based on the historical data. We demonstrate that the voluntary vaccination alone will not be enough to substantially limit the spread of MPX.

2 Model and methods

We adapt a SVEIR compartmental model of MPX [21] and extend it by incorporating the game-theoretic framework of voluntary vaccination as done in [71].

Individuals are born susceptible (S) at rate Λ. Without vaccination, a susceptible individual becomes exposed (E) after coming in contact with an infectious individual (I); this happens at rate βIN where N is the population size. The incubation period lasts σ^−1^ after which the individual develops MPX. The MPX cases recover (R) at rate γ^−1^ and gain permanent immunity.

The original model presented in [21] allowed for a proportion of individuals to be vaccinated at birth. However, since the smallpox vaccine is no longer mandated at birth, it is more realistic to consider that the vaccination occurs later in life. As in [65] for the COVID-19 vaccine, we will assume the susceptible individuals are vaccinated (V) at rate ψ. While in theory, ψ ∈ [0, ∞), there are bounds on how fast the population can be vaccinated. We will thus assume ψ ∈ [0, ψ~~max~~] where ψ~~max~~ is the maximal feasible vaccination rate.

The vaccine does not provide complete protection; the vaccine efficacy is e ∈ (0, 1]. It follows that the vaccinated individuals become exposed at rate (1-e)βIN.

Finally, all individuals can die of natural causes at a rate μ.

This compartmental model is illustrated in Fig 1. We note that this ODE model is a special case of a SARS model considered in [74]. The notation is explained and the parameter values are shown in Table 1. The model calibration is explained in detail in Section 4.

Fig 1: Scheme of the SVEIR model for MPX.

We extend this ODE model of MPX transmission by incorporating the game-theoretic component as done in [71] for measles, smallpox, and other childhood diseases or in [65] for the recent COVID-19 outbreak or [63] for an Ebola outbreak.

A vaccination game is played by susceptible individuals who are assumed to be rational, acting in their own best interests, and having complete information about the MPX epidemic. The individuals decide whether to vaccinate or stay unvaccinated. The payoff to the individual is a function that depends on the action of that individual (whether they vaccinate or not) and the actions of other players (how fast they are vaccinating as a whole). The payoff incorporates the cost of the vaccination, CV, the risk of getting infected, πNV and πV evaluated below, and the costs of the MPX infection CMPX*.

To evaluate the probability of getting exposed to MPX, we follow [71]. We assume that the epidemics reached a steady state with I* infected individuals. The formula for I* is given later by (15); it is important to note that I* depends on ψ, the vaccination rate in the population, but not on the decision of the focal individual.

The probability that an unvaccinated individual becomes exposed to MPX is πNV=βINβIN+μ, (1) where βIN+μ is a rate at which individuals with no intention to vaccinate leave the Susceptible compartment; βIN is the rate at which they enter the Exposed compartment.

Similarly, the probability that a vaccinated individual becomes exposed to MPX is πV=(1-e)βIN(1-e)βIN+μ. (2) Once exposed, the individual will become infected with probability σσ+μ.

The solution of the vaccination game, the Nash equilibrium (NE), is the population vaccination rate ψNE such that in this situation no individual has an incentive to deviate from the population strategy, i.e. either (1) ψNE = 0 when πNV < πV for ψ = 0, (2) ψNE = ψ if πNV = πV when the vaccination rate is ψ, or (3) ψNE = ψ~~max~~ if πNV > πV when the vaccination rate is ψ = ψ~~max~~.

3 Analysis

3.1 Analysis of the ODE system

The model yields the following differential dSdt=Λ-(ψ+μ+βIN)S (3) dVdt=ψS-(μ+(1-e)βIN)V (4) dEdt=(βNS+(1-e)βNV)I-(μ+σ)E (5) dIdt=σE-(μ+γ)I (6) dRdt=γI-μR. (7)

There are two equilibria of the dynamics (3)–(7). The disease-free equilibrium (DFE) E0=(S0,V0,E0,I0,R0) is given by (Λμ+ψ,Λψμ(μ+ψ),0,0,0), i.e., S0=Nμμ+ψ and V0=Nψμ+ψ where N=Λμ is the total population size. [74] derived the effective reproduction number, i.e., the number of secondary infections from a single infected individual in an otherwise healthy population, as R(ψ)=R0μ+(1-e)ψψ+μ. (8) Here, R0=σβ(σ+μ)(γ+μ) (9) is the basic reproduction number, i.e., the number of secondary infections from a single infected individual in an otherwise healthy and unvaccinated population.

The DFE is globally asymptotically stable if R(ψ) ≤ 1 [74, Theorem 4.1].

Let ψHI be the minimal level of vaccination needed for achieving a herd immunity; specifically let ψHI ∈ [0, ∞) be such that R(ψ)≤1 for all ψ ≥ ψHI. It follows that DFE is the only stable equilibrium for ψ ≥ ψHI and ψHI={0ifR0≤1,∞ife≤1-1R0,R0-11-(1-e)R0μotherwise. (10) In particular, if e≤1-1R0, then no vaccination rate will prevent the epidemic. Moreover, if ψHI > ψ~~max~~, then no feasible vaccination rate can prevent the epidemic.

The endemic equilibrium exists only when R(ψ) > 1. Moreover, the endemic equilibrium is locally asymptotically stable whenever R(ψ) > 1. The theoretical results and simulations performed in [74] suggest that it is globally asymptotically stable when R(ψ) > 1 and the initial population satisfies E|~~t = 0~~ > 0 or I|~~t = 0~~ > 0. These results are also supported by [75] who studied a similar model without the E compartment.

Thus, for any set of parameter values, there is only one stable equilibrium of the ODE system. Let us denote it by E*=(S*,E*,I*,R*,V*). As done in [74], setting the derivatives of (3)–(7) to 0 and solving the resulting system of algebraic equation yields a polynomial equation for I* in the form I*(aI2+bI+c)=0 (11) where a=(βN)2(1-e) (12) b=βNμ(1+(1-e)(ψ+μ)μ+(1-e)ψ(1-R(ψ)+(1-e)ψμ)) (13) c=μ(ψ+μ)(1-R(ψ)). (14) The equation aI*^2^ + bI* + c = 0 has no positive root when R(ψ) < 1. Thus, if R(ψ) ≤ 1, I* = 0 is the only biologically relevant solution of (11). On the other hand, as seen above, when R(ψ) > 1, the disease-free equilibrium corresponding to I* = 0 is not stable. Thus, we have I*={0,ifR(ψ)≤1,i.e.,ifψ≥ψHI,-cb,ifR(ψ)>1,ande=1,-b+b2-4ac2a,otherwise, (15) Note that I* is always non-increasing in ψ and it is decreasing in ψ whenever I* > 0 (and e > 0).

Furthermore, it follows easily from the algebra that S*=ΛβIN+ψ+μ (16) E=μ+γσI* (17) R*=γμI* (18) V*=ψS*(1-e)βI*N+μ. (19)

3.2 Solving the vaccination game

Let CMPX*=CMPXCV be the expected cost of MPX infection expressed relative to the cost of the vaccine. The incentive function h(ψ) for an individual to vaccinate when the population vaccination rate is ψ is given by a difference between the expected costs if not vaccinated and the expected costs if vaccinated (including the potential cost of contracting MPX), i.e., we can assume h(ψ)=CMPX*(σσ+μ)(πNV-πV)-1. (20) Here, (σσ+μ) is the probability an exposed individual becomes infected. It follows that the Nash equilibrium is given by ψNE={0,ifh(0)<0,root(s)ofh,ifin[0,ψmax],ψmax,ifh(ψmax)>0. (21)

3.3 Calculations of the Nash equilibria

Here we show detailed calculations for the Nash equilibria of the vaccination game and find the roots of h(ψ) = 0. We will study the function h˜(x)=(CMPXσσ+μ)(xx+1-(1-e)x(1-e)x+1)-1 (22) =(CMPXσσ+μ)xe(x+1)((1-e)x+1)-1. (23) Note that h(ψ)=h˜(βΛI*); and while we are primarily interested in the behavior of h(ψ) on [0, ψHI), we will investigate the function h˜(x) on [0, ∞).

It follows that h˜(x)⋚0 if and only if 0⋚(1-e)x2+[2-e(1+CMPXσσ+μ)]x+1. (24) We will assume CMPXσσ+μ>1, as otherwise the right-hand side of (24) is positive for x ≥ 0. Note that this is a reasonable assumption; μ ≪ σ and thus σ+μσ≈1 and if CMPX*<1, then the vaccine would cost more than the disease, i.e., nobody would vaccinate. We will also assume (2-e(1+CMPXσσ+μ))2-4(1-e)>0, i.e., 4CMPXσσ+μ(1+CMPXσσ+μ)2≤e<1. (25) When (25) does not hold, there are no real roots of h˜ and the right-hand side of (24) is positive for x ≥ 0. When (25) holds, the roots of h˜ are given by x1,2=e(1+CMPXσσ+μ)-2±(e(1+CMPX*σσ+μ)-2)2-4(1-e)2(1-e). (26)

Setting (15) equal to I1,2*=Λβx1,2 where x~~1,2~~ is given by (26), we obtain -μ2(1-e)x1,22=μ(ψ+μ)(1-R(ψ))+μ2(1+(1-e)(ψ+μ)μ+(1-e)ψ(1-R(ψ)+(1-e)ψμ))x1,2. (27) Simplifying (27) and solving for ψ yields ψ1,2=μ(-(1-e)x1,22+((1-e)(R0-1)-1)x1,2+R0-1)1+(1-e)(x1,2-R0). (28)

3.4 Analysis of the Nash equilibria

In the context of vaccination games studied in this paper, there are two kinds of NE. If h(ψNE) = 0 and h′(ψNE) < 0, then ψNE is called a convergent stable Nash equilibrium (CSNE) [71]. If the population adopts a strategy ψ ≈ ψNE, then it will evolve even closer to ψNE. For similar reasons, if h(0) < 0 or if h(ψ~~max~~) > 0, then 0 or ψ~~max~~ are CSNEs. However, when h(ψNE) = 0 and h′(ψNE) > 0, then ψNE is not CSNE, as a small deviation from ψNE will result in even larger deviation.

We note that there is a difference between the usual vaccination games (such as [70, 71]) that assume e = 1 and our more general case with e < 1. If e = 1, the incentive function h(ψ) is decreasing in ψ and ψHI < ∞. Thus, assuming ψHI < ψ~~max~~, there is ψ > ψHI for which I* = 0 and h(ψ) = −1 < 0. Consequently, there is at most one root of h(ψ) which exists if and only if h(0) ≥ 0. However, when e < 1, there can be multiple roots of h and also multiple Nash equilibria as demonstrated later.

3.5 Uncertainty and sensitivity analysis

We performed uncertainty and sensitivity analysis using the Latin hyper-cube sampling with partial rank correlation coefficient (LHS-PRCC) scheme [76, 77]. The scheme is described in detail in [78] and their MATLAB and R implementation can be found in [79]. Our MATLAB code, including the code for uncertainty and sensitivity analysis, is in the S1 Code.

The Latin Hyper-cube Sampling (LHS) is a stratified sampling without replacement; the random parameter distributions are divided into intervals of equal probability and the sampling is independent for each parameter. This method provides an unbiased estimate of the average model output while it requires fewer samples than simple random sampling to achieve the same accuracy [80].

The partial rank correlation coefficient (PRCC) between model parameter p and model output O is a correlation coefficient rRp,RO=Cov(Rp,RO)Var(Rp)Var(RO) (29) between Rp and RO which are residuals of the rank-transformed linear regression models for p and O. PRCC is a robust sensitivity measure for nonlinear but monotonic relationships between inputs and the output, as long as little to no correlation exists between the inputs [78].

4 MPX calibration

4.1 Demographic parameters

We used data from the CDC to establish the natural birth rate and natural death rate. From [73], the birth rate in the U.S. in 2020 was 11.0 per 1,000. We will thus assume μ=111000·365 per individual per day. From [72], the death rate in the U.S. in 2020 was 1,027.0 deaths per 100,000 population. We will thus assume Λ=0.01027365.

4.2 Disease progression

The MPX incubation period lasts 7–17 days [2]. We thus assume σ^−1^ = 12 days on average.

The prodromal period lasts 1–4 days [2] and the rash period lasts 14–28 days [2]. We thus assume that γ^−1^ = 23.5 days on average.

4.3 MPX related mortality

For the simplicity of the ODE model and its analysis, we do not consider any MPX related mortality. While the mortality was reported to be as high as 11% [4], and in recent times the case fatality has been 3–6% [81], the current 2022 outbreak has over 6000 reported cases of MPX and only 3 deaths [6]. We note that [82] reports 66 deaths in African countries and there may be a time lag in death reporting. Thus, there may be a non-negligible mortality even in the current outbreak. Yet, the mortality is likely relatively small to make substantial impacts on the main conclusions.

4.4 Vaccine efficacy

Data from Equateur Province of Democratic Republic of Congo from 1981 to 1986, in the years following smallpox eradication, suggest that smallpox vaccine conferred 85% protection against monkeypox [2]. For the purpose of this paper, we will thus assume e ≈ 0.85.

The data are described in detail in [4]. Within the household contacts, the attack rate among unvaccinated individuals was 9.28 and among the vaccinated individuals it was 1.31 [4, Table 5]. Thus, the vaccine efficiency was estimated as 1-1.319.28≈0.855. However, the numbers within the house differ by gender. For males, the efficacy was 1-0.858.61≈0.9 while for females the efficacy was 1-1.749.91≈0.82.

Given that in the 2022 MPX outbreak, 99.5% (4385/4406) of cases for which the sex is known are men [6], it may be tempting to adopt e = 0.9 based on the above calculations. However, the nature of contacts for data collected in [4] was likely different from the nature of contact in the 2022 outbreak in MSM.

The 3rd generation smallpox vaccine Imvanex (Modified Vaccine Ankara—MVA) has been authorised by the European Medicines Agency, but scientific evidence on the vaccine effectiveness of MVA against MPX is still lacking [84]. Overall, a more recent estimate of vaccine efficacy is needed.

4.5 Basic reproduction number

While the basic reproduction number is not a parameter of the model, we will use the estimate of R0 to derive an estimate of β by R0γ as in a section outlined below.

Using data collected in the Democratic Republic of the Congo between 1966–1984, [21] estimated the basic reproduction number for monkeypox as R0 ≈ 2.13 with bounds between 1.46 and 2.67.

To estimate R0 for the 2022 MPX outbreak, we used data from [85], who shared their raw data set of MPX cases on github [86]. We used SAS to obtain incidence rates for the whole world as well as for Spain, England, and Germany, which were the three countries with the most cases; the SAS code is provided in the S2 Code.

We used an online tool [83] to estimate R0 from the incidence data as follows. We set the sliding window to 7 days, set prior mean value for R0 to 2.1, and prior SD to 0.3. We assumed all transmissions to be local (for the lack of better data). We used the serial interval by distributional estimate, using the option “Parametric without uncertainty” and set the mean to 9.7, and SD to 0.5 [87]. The results are summarized in Table 2. While the estimates vary and generally decrease in time, there are some relatively high values of R. This leaves the possibility of R0 > 3 open at least for the current 2022 outbreak.

4.6 Transmission rate

By (9), and assuming μ ≪ σ and μ ≪ γ, R0=σβ(σ+μ)(γ+μ)≈βγ. (30) Thus, β≈R0γ. (31) Based on estimates of R0 from [21] discussed above, we get β ≈ 2.13/23.5 ≈ 0.09 with bounds between 1.47/32 ≈ 0.045 and 2.67/14 ≈ 0.18. We note that this estimate agrees with a relatively crude estimate from [70] who used transmission risk data from [88] to arrive at the same value.

4.7 Validation

To validate our choice of parameters and their ranges, we run the sensitivity analysis of the basic reproduction number on the parameter values and ranges as specified above and in Table 1.

The estimated basic reproduction number is R0 ≈ 2.11, and the average value during the uncertainty analysis was 2.53; see Fig 2. This value seemed reasonable and in agreement with historical data; see Fig 3.

Fig 2: Uncertainty analysis of R0.Parameter values and ranges as specified in Table 1. The average value of R0 is about 2.53.

Fig 3: Actual (black dots) and predicted incidence.Full red line is for β = 0.09, dashed is for β = 0.19 and dotted is for β = 0.29. All other parameters are as specified in Table 1.

We also plotted the actual new MPX cases as obtained from [85] against the model’s predicted incidence. We numerically solved the system (3)–(7). The number of new cases at day d was obtained as σE|~~t=d~~. We normalized it so that at the start of the epidemic we have σE|~~t=0~~ = 1; see Fig 3. We note that the best match appears for β ≈ 0.29. However, this would yield R0 ≈ 6.8 which is more in the order of R0 for smallpox [89]. The fit is still good enough for β ≈ 0.18 which yields R0 ≈ 4.23. This is still larger than the usual estimates for MPX. However, as shown in Table 2, such a value of R0 is not completely unreasonable for the 2022 MPX outbreak. Consequently, while we will do most of our calculations for β = 0.09 which seems to be in agreement with all previous estimates and historical data, we will also consider β = 0.18, the upper bound for β estimates.

4.8 Maximal feasible vaccination rate

We did not locate any data on the maximal feasible vaccination rate. We will assume that the population can get vaccinated about once a year, i.e., ψ~~max~~ = 1/365. In the U.S., at the time of writing, the demand for the vaccine exceeds the supply [90] and for large populations, the vaccine supply is likely the most important factor limiting the vaccination rate. However, even for smaller populations and/or during the times the vaccine supply will be restored, there may be logistical issues (such as limited supply of qualified nurses and doctors) preventing a significantly faster vaccination rate.

4.9 Costs of vaccination and costs of MPX

There can be many types of costs associated with vaccination, including the actual cost for vaccination, time loss, and travel cost that all have negative effects on the probability of complete vaccination [91].

Similarly, the expected cost of the disease C~~MPX~~ includes possible medication costs, doctor charges, time loss and similar direct and indirect costs.

We were not able to locate any reliable and accurate values for CV and C~~MPX~~ applicable to the current outbreak. We note that [70] used CV = 4 based on [92] and C~~MPX~~ = 100 based on [93]. However, those values are for the Democratic Republic of the Congo and based on historical data. In the current 2022 outbreak, especially in the U.S., the vaccine is provided for free, and the cost is thus limited to indirect costs such as taking the time off work to get vaccinated and possibly dealing with minor vaccine side effects [94]. At the same time, the MPX infection does not seem to add any extra monetary expenses to the individuals, apart from mild symptoms and taking the time off from work although cases of severe pain have been reported. We will thus assume CMPX* to be somewhere between 1 and 10 and note that in reality, the value will likely be different from person to person.

The value of CMPX* plays a role in the incentive function given in (20). The formula contains (σσ+μ), the probability an exposed individual becomes infected, which, for MPX parameters discussed above, is approximately 1. Here, we can assume that CMPX* is not affected by natural mortality. However, if the MPX infection lasts significantly longer, one would have to carefully account for the effects of natural mortality.

Furthermore, we note that while the real costs are important, what truly matters for the individuals are the perceived costs of the disease and vaccination; yet the actual model does not change when we change the interpretation of CMPX* from “real” to “perceived” relative cost.

5 Results

Historical data indicate that MPX outbreaks can be eliminated by vaccination. Indeed, the estimates for the vaccine efficacy e ≈ 0.85 and basic reproduction number R0 ≈ 2.13 based on data from 1966–1984 [21] mean that e>1-1R0. Thus, the vaccination rate needed to achieve herd immunity, ψHI, is finite. For the parameter values considered in Table 1, we get μ/ψHI ≈ 0.61, meaning that if the whole population at risk can be vaccinated in about 61% of the average lifespan, then herd immunity will be reached. Fig 4 illustrates the uncertainty and Fig 5 shows the sensitivity analysis of μ/ψHI. It shows a natural result that in order to achieve MPX elimination, the vaccination has to be done faster if the transmission rate, β, or the infectious period, γ^−1^, increase. On the other hand, the vaccination can be done slower if the vaccine efficacy, e, increases.

Fig 4: Uncertainty analysis of μ/ψHI.Parameter values and ranges as specified in Table 1. For those values, μ/ψHI ≈ 0.61 meaning that to achieve herd immunity, the whole population at risk has to be vaccinated in about 61% of the average lifespan. The average value of μ/ψHI is about 0.53.

Fig 5: Sensitivity analysis of μ/ψHI.Parameter values and ranges as specified in Table 1. For those values, μ/ψHI ≈ 0.61 meaning that to achieve herd immunity, the whole population at risk has to be vaccinated in about 61% of the average lifespan.

Figs 6–8 illustrate the Nash equilibrium values, the MPX prevalence, and the annual MPX incidence in the population using the optimal voluntary vaccination rates. For the parameters as in Table 1, the NE is to “not vaccinate” as long as CMPX*, the cost of MPX relative to the cost of vaccine, is less than about 2.6. Without vaccination, the prevalence would be about 3.5 cases per 10^4^ individuals. However, even with optimal voluntary vaccination and relatively high CMPX*≈10, the MPX prevalence in the equilibrium is more than 0.5 cases per 10^4^ population. The annual incidence without any vaccination is almost 60 cases per 10^4^; and even with optimal voluntary vaccination and CMPX*≈10, the incidence would stay around 10 cases per 10^4^ population.

Fig 6: The optimal voluntary vaccination rates.The lines are color coded corresponding to the regions shown in Fig 9 which shows a diagram for Nash equilibria as e and CMPX* vary. Blue: 0 is the only NE and it is CSNE. Brown: positive ψNE is the only NE and it is CSNE. Red dashed line shows the value of ψHI. Unless varied or otherwise specified, the parameters are as in Table 1.

Fig 7: MPX prevalence in a population that uses optimal voluntary vaccination rates.The lines are color coded corresponding to the regions shown in Fig 9 which shows a diagram for Nash equilibria as e and CMPX* vary. Blue: 0 is the only NE and it is CSNE. Brown: positive ψNE is the only NE and it is CSNE. Unless varied or otherwise specified, the parameters are as in Table 1.

Fig 8: Annual MPX incidence in a population that uses optimal voluntary vaccination rates.The lines are color coded corresponding to the regions shown in Fig 9 which shows a diagram for Nash equilibria as e and CMPX* vary. Blue: 0 is the only NE and it is CSNE. Brown: positive ψNE is the only NE and it is CSNE. Unless varied or otherwise specified, the parameters are as in Table 1.

To better illustrate how the outcomes depend on different parameter values, Fig 9 shows the NE as e and CMPX* vary while β = 0.09 (and thus R0 ≈ 2.11). For these values, there are only two possibilities, either ψNE = 0 is the only NE or there is a unique positive NE, 0 < ψNE < ψ~~max~~. Fig 10 further illustrates what happens when the transmission rate β increases to 0.18, or equivalently, if R0 increases to 4.23. There are now four distinct the region where 0 is the only NE is almost the same, but the region with unique positive ψNE < ψ~~max~~ is significantly smaller and a new region where ψNE = ψ~~max~~ appeared for medium values of e and large enough values of CMPX*. More importantly, there is also a region with three distinct NE. Fig 11 shows the regions of NE when CMPX* and R0 varies while e = 0.85. We can see that the region with three NE exists only for relatively high values of R0 (about 4 and more). The region exists only for a narrow range of values of CMPX*, although the range gets wider as R0 increases.

Fig 9: Nash equilibria as e and CMPX* vary.Blue: 0 is the only NE and it is CSNE. Brown: positive ψNE is the only NE and it is CSNE. Unless varied or otherwise specified, the parameters are as in Table 1.

Fig 10: Nash equilibria as e and CMPX* vary for higher transmission rate, β = 0.18.Other parameters are as in Table 1 unless they vary or are otherwise specified. Blue: 0 is the only NE and it is CSNE. Brown: positive ψNE < ψmax is the only NE and it is CSNE. Light three NEs, 0 and the larger NE are CSNE. Red: maximal feasible vaccination rate is the only CSNE.

Fig 11: Nash equilibria as R0 and CMPX* vary; β is estimated by (21) as β ≈ R0γ.Other parameters are as in Table 1. Blue: 0 is the only NE and it is CSNE. Brown: positive ψNE < ψmax is the only NE and it is CSNE. Light three NEs, 0 and the larger NE are CSNE. Red: maximal feasible vaccination rate is the only CSNE.

Figs 12–15 show graphs of incentive functions corresponding to the different NE regions in Figs 9–11. Specifically, Fig 12 illustrates an incentive function in a region where ψNE = 0 is the only NE. Fig 13 is an example of an incentive function where ψNE ∈ (0, ψ~~max~~) is the only NE. The incentive function with multiple roots is shown in Fig 14, while the incentive function for the region where NE is the maximal vaccination rate is shown in Fig 15.

Fig 12: The incentive function for the parameters in the blue region of Fig 9 where 0 is the only NE and CSNE.Full black circle is the CSNE, the red circle corresponds to ψHI. Unless varied or otherwise specified, the parameters are as in Table 1.

Fig 13: The incentive function for the parameters in the brown region of Fig 9 where ψNE > 0 is the only NE and CSNE.Full black circle is the CSNE, the red circle corresponds to ψHI. Unless varied or otherwise specified, the parameters are as in Table 1.

Fig 14: The incentive function for the parameters in the light blue region of Fig 10.There are three NE at the same time. Full black circles are the CSNE. The empty circle is NE that is not CSNE. The red circle corresponds to ψHI. Unless varied or otherwise specified, the parameters are as in Table 1.

Fig 15: The incentive function for the parameters in the red region of Fig 10.Full black circle is the CSNE. There is no ψHI. Unless varied or otherwise specified, the parameters are as in Table 1.

Figs 16–18 illustrate what happens when the transmission rate β increases to 0.18. Without vaccination, the MPX prevalence and annual incidence increase roughly by a factor of 1.5. There is also a backward bifurcation for CMPX*≈2.25. For those values, we have three optimal voluntary vaccination rates. Only 0 and the largest value are CSNE. The medium vaccination rate, ψ~~NE,1~~ is NE but not CSNE. From the public health perspective, it means that there is a need for public policy to mandate the vaccination rate to be at least ψ~~NE,1~~; otherwise the voluntary rate would decline to 0. When CMPX* is large enough for ψNE > 0 to exist, the vaccination rate is larger in this case compared to case when β = 0.09. Consequently, the prevalence and incidence in the population using ψNE is lower when β is bigger.

To further assess what could happen without vaccination, Figs 19 and 20 show the uncertainty and sensitivity analysis of MPX prevalence and incidence of unvaccinated population. The average prevalence is around 6.5 cases in 10^4^ population and the average annual incidence is around 110 cases per 10^4^ population, further underlying the importance of vaccinations to try to curtail the outbreak.

Fig 19: Uncertainty analysis for the MPX prevalence in 10^4^ unvaccinated population (the average is approximately 6.5).

Fig 20: Uncertainty analysis for the MPX annual incidence per 10^4^ unvaccinated population (the average is approximately 110).

6 Conclusions and discussion

We applied the vaccination game theory framework developed by [71] to the compartmental model of MPX transmission [21] explicitly incorporating the possibility of MPX infections even for the vaccinated population.

Without vaccination, MPX could become endemic with relatively high prevalence (3.5 cases per 10^4^) and incidence (almost 60 cases per year per 10^4^) levels. We identified optimal voluntary vaccinations rates, i.e. rates that are likely to be adopted by the population without any central or government mandates and interventions. For a relatively low cost of MPX infection (less than 2.5 times the cost of the vaccine), to not vaccinate is unfortunately an optimal strategy from the individual standpoint. Even as the cost of infection increases, the optimal voluntary vaccination rate is not enough to substantially decrease the number of MPX cases. Mandatory vaccination for individuals at risk is therefore highly recommended.

The result that voluntary vaccination alone is not enough to eliminate MPX is not surprising. It is a consequence of the low cost of MPX infection (relative to the cost of vaccination). It has already been demonstrated before that the tendency of individuals to optimize self-interest can lead to vaccination levels that are suboptimal for a community [95] and similar predictions have been made in general [96] as well as for specific scenarios involving yellow fever [97], typhoid fever [98], cholera [99], and Hepatitis B [100].

Moreover, we demonstrated that, for some parameter values, specifically a relatively low cost of MPX infection and relatively high rate of MPX transmission, there are multiple Nash equilibria of the vaccination game and the solutions exhibit backward bifurcation. For the public health officials, this means that a minimal vaccination rate has to be mandated in this case, as otherwise the population vaccination rate would decline to 0. The existence of multiple Nash equilibria for vaccination games is a relatively new and not yet fully investigated phenomenon. To our knowledge, [65] is the only other work in the vaccination game theory where multiple equilibria occur for a single action, although the backward bifurcation has not been investigated there; moreover [97] and [101] investigated multiple equilibria in vaccination games with two preventive actions.

Our results underline the importance of proper estimation of the vaccine efficacy and the reproduction number for the current MPX outbreak [102]. The estimates of basic reproduction number R0 ≈ 2.13 [21], the effective reproduction number R ≈ 0.83 [8, 103] and vaccine efficacy e ≈ 0.85 [8] are all based on historical data from 1966—1984, during or soon after smallpox vaccinations ended. The current outbreak affects primarily MSM, not children, and there are signs that R0 can be bigger than expected [104]. If R0 is larger than 4, then even a slight decrease of vaccine efficacy can mean that even with full vaccination, the MPX outbreak may not be stopped. Moreover, in that scenario, the transmission rate would likely be large enough to have multiple Nash equilibria and backward bifurcation.

As with any other mathematical model, our model has a number of limitations and simplifying assumptions.

We performed the analysis as if MPX already reached the endemic state, which is fortunately not yet the case for most of the countries. The vaccination adoption behaviour can happen at about the same time scale as the infection dynamics, allowing for co-evolution [105]. The coupling of game and epidemic models can lead to oscillations in vaccine uptake over time [95]. The vaccine-generated herd immunity can lower disease incidence so much that real or perceived vaccine risks causes individuals to cease vaccinating which in turn causes uptick in disease incidence [106]. There is even a potential for a significant instability if the perceptions of vaccine and infection risks are homogeneous in the population [95]. We expect that the oscillations would be even more profound with the vaccine imperfection.

Furthermore, we assumed that individuals are well informed about MPX which is also not the case. In the U.S., almost half the respondents (47%) feel that their knowledge level about MPX is poor or very poor [107]. Moreover, infectious diseases can be under-reported, and MPX is no exception with testing only recently expanding [108]. Also, at the time of writing, the demand for the vaccine exceeds the supply [90], i.e. the maximal feasible vaccination rate is relatively low. There is now a growing body of literature on disease transmission and misinformation [109–111]. Misinformation can prevent the suppression of epidemics [112]. It is quite conceivable that, similarly to what happened with HPV vaccine in Denmark [113], a misinformation about MPX vaccine can diminish vaccine coverage.

Perhaps the most severe limitation is that we assumed homogeneous well mixed population and, as a result, we obtained a single Nash equilibrium for most parameter values. Complex networks provide a better platform for more realistic modeling [52, 114–116] and explicitly incorporating social networks within MSM community [117–119] would thus greatly improve the model. Heterogeneity in the population yields the heterogeneity in vaccinating actions [120]. The individuals with many contacts typically have higher inclination to voluntary vaccinate and this can help inhibit the outbreaks [121]. As another extension into heterogeneous populations, one could incorporate the fact that different individuals can perceive the cost of MPX and the cost of vaccination differently. This assumption could significantly alter our results. In the present (homogeneous) model, if the cost of MPX is smaller than the cost of vaccination, individuals opt not to vaccinate. However, in the heterogeneous model, even if the average cost of MPX is smaller than the average cost of vaccination, there can still be a non-negligible proportion of population who perceive the cost of MPX as significantly larger than the cost of vaccine and, as a result, opt to vaccinate.

One has to keep in mind that no model can be fully realistic and account for every detail [122, 123] and so, despite all the above shortcomings, our model provides a reasonable, and not so positive, outlook into what could happen without any mandate for vaccinations and/or possibly other measures to stop the unfolding MPX outbreak.

Voluntary vaccination may not stop monkeypox outbreak: A game-theoretic model