Vaccination games in prevention of infectious diseases with application to COVID-19.

Vaccination coverage is crucial for disease prevention and control. An appropriate combination of compulsory vaccination with voluntary vaccination is necessary to achieve the goal of herd immunity for some epidemic diseases such as measles and COVID-19. A mathematical model is proposed that incorporates both compulsory vaccination and voluntary vaccination, where a decision of voluntary vaccination is made on the basis of game evaluation by comparing the expected returns of different strategies. It is shown that the threshold of disease invasion is determined by the reproduction numbers, and an over-response in magnitude or information interval in the dynamic games could induce periodic oscillations from the Hopf bifurcation. The theoretical results are applied to COVID-19 to find out the strategies for protective immune barrier against virus variants.

Introduction

The spread of infectious diseases, for example, highly infectious measles, hepatitis B and the recent epidemic of COVID-19, poses a great threat to our lives and to economic development in all regions of the globe [1], [2], [3]. Though the prevention measures such as early detection, early quarantine and early treatment are important, vaccination is the most effective to control infectious diseases. In fact, the highly contagious smallpox and polio were eradicated or eliminated through vaccination campaigns [4], [5], [6]. The success of a vaccination program is crucially dependent upon the vaccination coverage above which, the spread of infectious disease is controlled. The approach to fulfill that goal is affected by multiple factors, including the attitude of individuals for vaccination, infection risk of disease and side-effect of vaccine. Therefore, it is necessary to find an optimal way to deploy the vaccination program.

An appropriate fraction of compulsory vaccination for individuals with the higher infectious transmissions is important to curb the spread of infectious diseases. For example, employees in public services such as health-care workers, bus drivers and teachers have the much more probability to be infected and to transmit the disease. Therefore, the compulsory vaccination for the leading infectious group is fundamental in a vaccination program. However, the mandatory vaccination is limited by the beliefs of free choices in vaccination or religious faith against vaccination. Thus, the voluntary vaccination applies to the majority of a population. This means that the coverage of voluntary vaccination is critical in order to achieve the goal of immune barrier. However, for voluntary vaccination, the decision is made on the basis of evaluating the cost and benefits of vaccination. For example, the decision is given by comparing the loss from a potential infection with the safety and cost of vaccination [7], [8], [9].

Game theory is often applied to vaccination decisions, in which individuals choose the most advantageous strategy by comparing the expected returns of two strategies. The researches focus upon the case where all the individuals in a population are voluntarily vaccinated and they respond linearly to the pay-off information [10], [11], [12], [13]. Moreover, since a vaccination decision may be delayed at a certain time, a discrete time delay or a compartment for the delayed vaccination are incorporated into the dynamic game models [14], [15]. Notice that the accumulated historical information of disease status plays the key role in a realistic vaccination decision [16], [17]. Furthermore, a nonlinear response function for vaccination decision is appropriate in vaccination games [18], [19], [20]. In this paper, we propose a mathematical model with the vaccination game, which incorporates the integral information of disease status and introduces a nonlinear response function in decision-making game. Another novelty of the model is that both mandatory vaccination and voluntary vaccination are included, which can be used to analyze the impact of vaccination behaviors on the prevention and control of COVID-19.

The organization of this paper is as follows. In the next section, we formulate the mathematical model of epidemic diseases with behavior changes for vaccination. Section 3 presents the mathematical analysis for the stability and Hopf bifurcation of equilibrium points of the model. In Section 4, we use the model to simulate the dynamics of COVID-19 spread in China and analyze the vaccination strategies for immune barrier. The paper ends with some discussions in Section 5.

Model

We start with the formulation of mathematical model of epidemic diseases of SIRVS type, which is applicable to describe the transmission dynamics of infectious diseases such as measles, influenza, haemorrhagic fever and COVID-19. Here, the population is simply divided into four groups: susceptible (S), vaccinated (V), acutely infected (I), and recovered (R) individuals. The individuals with vaccination failure become susceptible, susceptible individuals enter the infection compartment after valid infectious contacts with infectors, and infected individuals are admitted to the recovery compartment due to natural recovery, quarantine and treatment. Notice that reinfections in vaccinated individuals or naturally infected individuals are common in the diseases such as COVID-19 and influenza [21], [22], [23]. We assume that recovered individuals can return into the susceptible compartment. The flowchart of disease transmission and progression is shown in Fig. 1 , and the dynamics of the state variables are described by the following differential equations:where θ is the proportion of susceptible people who are compulsorily vaccinated because they work in high-risk occupations, such as doctors, teachers and staffs in transportation services; other individuals are voluntarily vaccinated with the vaccination proportion x  (0 ≤ x ≤ 1); σ represents the vaccine efficacy (0 ≤ σ ≤ 1); ϕ is the waning rate of immunity in vaccinated individuals; κ is the waning rate of immunity in recovered individuals; β is the effective transmission rate; γ is the natural recovery rate. Furthermore, q = q 0 ⋅ γ 2 where q 0 is the quarantine rate and γ 2 is the treatment rate after quarantine. In addition, we use the coefficient ρ  (0 < ρ < 1) to describe the limited supply of vaccine.

Fig. 1

The transmission and progression of infectious disease.

Note that the birth rate and death rate in the population are identical and the disease-induced death is ignored. It is easy to see that the population size is a constant N. Let s = S/N, i = I/N, r = R/N, v = V/N so that s + i + r + v = 1. Then (2.1) is simplified towhere the equation in r is decoupled.

Notice that the vaccination proportion x is altered by vaccination behaviors, which is a dynamical decision-making game. To describe such a game, we assume that the cost of vaccination (including the cost of vaccine, the burden of adverse reactions after vaccination, etc.) is c and the cost of infection (including the cost of treatment and the loss) is c . According to the dynamics of disease transmission, the rate at which an individual is infected is βi. Thus, the infection risk for an unvaccinated individual could be measured by βi ⋅ c . Moreover, (1 − σ)βi ⋅ c is the infection risk of an individual with unsuccessful vaccination, and σϕ ⋅ βi ⋅ c is the infection risk of an individual with temporary protection from vaccination. Let f be the pay-off for a vaccination strategy and let f be the pay-off for a non-vaccination strategy. Then we have

Let g 1(t) = f − f be the first information function that one defines the vaccination cost. The second one is the speed of new infections, i.e., the number of new cases per day, which is g 2(t) = βNsi. This function characterises the importance in change of disease prevalence. It is assumed that the net information function for decision-making is given bywhere m is a weight coefficient. If m = 1, the decision is solely dependent on the vaccination cost. For the case m = 0, the decision is made by considering only the speed of new infections.

It is assumed that individuals make vaccination decisions by considering both the information of present moment and the information of past time. In the aim to mimic the situation that the accumulated information over a time interval is crucial for vaccination decision, we take the exponential fading memory function K(t) = αe − and integrate the information to obtainwhere ε > 0 reflects an individual's sensitivity to current information and τ is the length of information interval.

Note that the value of decision-making function varies with disease status, which affects the inclination of individuals for vaccination. Motivated by previous studies [18], [24], we define the behavior response bywhere b > 0 represents the response intensity. This function is a continuous counterpart of classical Fermi function, which describes behavior switches in discrete imitation dynamics. Then following the classical approach to model imitation dynamics [12], [13], [16], [17], [19], we get the dynamical equation for behavior changes:where ν > 0 is a coefficient for imitation intensity.

Combining the disease transmission dynamics with Eq. (2.6), we get the full mathematical model:

Analysis

System (2.7) has a disease-free equilibrium with no voluntary vaccinationand a disease-free equilibrium with full voluntary vaccination

By the computation methods in papers [25], [26], we obtain the basic reproduction number of model (2.7):

The effective reproduction numbers at disease-free equilibrium points E 1 and E 2 are

Setwhere

If R 0 > 1, there is a boundary equilibrium point:

Furthermore, if R 1 > 1, there is an equilibriumwhere

All of the four equilibrium points are the boundary equilibrium points. Let us look for an internal equilibrium point (s 5,  i 5,  v 5,  x 5) of system (2.7). Setting the right sides of equations in system (2.7) to 0, we obtainand M = 0. The second equation of (3.3) implies s 5 = 1/R 0. From M = 0 we get g = 0, which implies

Then it is easy to obtain

Set

Evidently, (3.4), (3.5) mean

Assume

It follows from (3.5) that

Notice that

We have s 5 < 1. Using (3.3), we getwhich implies

Consequently, we can state the conclusion for the existence of interior equilibrium in (2.7).

Model(2.7)admits a unique interior equilibrium E5 = (s 5,  i 5,  v 5,  x 5) if (3.6) is satisfied.

We turn to consider the local stability of the equilibria. Set E = (s ,  i ,  v ,  x ) for n = 1,2,3,4,5 and

Linearize the system (2.7) at E and rewrite the state variables as s, i, v and x to obtainwhere

Let be an exponential solution of the linearized system, where λ and are constants, and substitute it into the linear system. Usingwe obtain the characteristic equation at E where

By analyzing the roots of the characteristic equation, we obtain the conditions for the local stability of the boundary equilibria, which are shown in Table 1 .

Table 1

The stability of boundary equilibrium points.

Equilibriums	Stability	Stability conditions
E1	Locally asymptotically stable	Rv0 < 1
E2	Always unstable
E3	Locally asymptotically stable	1 < Rv0 < Rc
E4	Locally asymptotically stable	Rv1 > Rc

For brevity, we present only the stability analysis for equilibrium E 3. First, λ 1 = P 3(3) is a characteristic root. The other three characteristic roots solve the equation:

Here,

Direct calculations indicate

It follows from the Routh-Hurwitz criterion that E 3 is asymptotically stable if

Note thatis equivalent to

Since the existence of E 3 implies R 0 > 1, we conclude that E 3 is asymptotically stable if

Let us now analyze the stability of equilibrium E 5. For the case of τ = 0, the characteristic equation of the linearized system at E 5 iswherein which,

By the Routh-Hurwitz criterion, we see that all the characteristic roots have negative real parts ifwhich implies that equilibrium E 5 is locally asymptotically stable when τ = 0. Motivated by numerical calculations in Example 3.1, we assume that there is a ν = ν 0 such that

Then there are a pair of pure imaginary roots λ 1, 2 = ± iω 0 in Eq. (3.9). In addition, from Eq. (3.9) one obtains

Hence, the other two roots of (3.9) have negative real parts. Note that

Consequently, we can state the following result.

If there is a ν = ν 0 such that Eq. (3.11) holds and . Then system (2.7) has a Hopf bifurcation at ν = ν 0 when τ = 0.

For the case of τ > 0, we consider below the existence of the Hopf bifurcation for the special case ε = 0. The analysis of the case ε > 0 is similar and will not be described here. The characteristic equation of the linearized system at E 5 is

Here,

In the aim to obtain the conditions under which E 5 keeps the stability for all τ > 0, we set λ = i ω in Eq. (3.12), and then separate the real and imaginary parts of the equation to get two transcendental equations:

Let z = ω 2, square both sides of the transcendental equation respectively and add the two transcendental equations to obtain:where

Set

Suppose that Eq. (3.14) has a positive root z. Then ±i ω with are a pair of pure imaginary root of characteristic equation. It follows from the transcendental equations that τ = τ (,  j = 1, 2, ⋯, are given by

Setwhere z ∗ is a root of Eq. (3.14) at τ = τ 0. By [27], [28], [29], all the roots of the characteristic eq. (3.12) have negative real parts on τ ∈ [0,  τ 0), which means that the equilibrium point E 5 is locally asymptotically stable.

In order to verify the cross-sectional condition, we differentiate both sides of Eq. (3.12) with respect to τ to getwhere

If h′(z ∗) ≠ 0, we have . Consequently, by [27], [28], [29], [30] we obtain the following conclusions.

Assume that there exist τ0 > 0 and ω 0 > 0 which satisfy (3.16) . Then equilibrium E 5 is locally asymptotically stable on [0,  τ 0). In addition, if h′(z ∗) ≠ 0, the system (2.7) has a Hopf bifurcation at τ = τ 0 .

We now illustrate the bifurcations of model (2.7) through numerical simulations.

Fix the parameters by β = 0.5,  σ = 0.9,  ϕ = 0.05,  μ = 0.05,   q = 0.01,  m = 0.85  ρ = 0.6,  γ = 0.1,  c = 10, c = 3000,  b = 8,  θ = 0.37,  α = 0.4,  κ = 0.001,  N = 1000. Then the model has an internal equilibrium point E 5 = (0.32,0.0076,0.6559,0.0152). First, we select τ = 0,  ε = 0.4. Then the equilibrium is asymptotically stable when ν < 0.0617. At ν = 0.0617, condition (3.11) is satisfied and . By Theorem (3.2), a Hopf bifurcation emerges at ν = 0.0617. Second, we choose ν = 0.02. If ε = 0, then τ = 4.6 and z ∗ = 0.0000436 satisfy Eq. (3.13) and Eq. (3.14), and h′(z ∗) = 1.85 × 10−8 ≠ 0. It follows from Theorem 3.3 that a Hopf bifurcation emerges at τ = 4.6. Finally, if ε = 0.2, a Hopf bifurcation emerges at τ = 5.7. These conclusions are verified by the plots of solutions of model (2.7) in Fig. 2 .

Fig. 2

Graphs of Hopf bifurcations for ODE case and DDE case in model (2.7).

Fix the parameters by μ = 0.1,  β = 0.6,  γ = 0.12,  ρ = 0.4,  ϕ = 0.001,  c = 20,  c = 1000,  m = 0.5, q = 0.001,  N = 1000,  θ = 0.55,  κ = 0.001,  b = 8,  ν = 2,  τ = 0,  ε = 0.4. Then E 2 is always unstable. The existence and stability of the other equilibria are shown in Fig. 3 . In particular, the system has five equilibrium points when 0.3644 < σ < 0.433, and has periodic solutions if 0.3644 < σ < 0.63.

Fig. 3

Bifurcation diagram of equilibria with the change of σ. The red, blue, pink and green curves correspond to the equilibrium points E1, E3, E4 and E5, respectively. In addition, the solid (dotted) lines indicate that the corresponding equilibrium points are stable (unstable), and H represents the Hopf bifurcation point.

Fix the parameters by μ = 0.1,  β = 0.75,  γ = 0.1,  ρ = 0.3,  ϕ = 0.005,  c = 2000,  m = 0.85,  q = 0.01, N = 1000,  θ = 0.6,  κ = 0.001,  b = 8,  ν = 0.07,  τ = 0,  ε = 0.4,  σ = 0.87. The parameter c does not affect the existence and the size of the boundary equilibrium points, but influences their stability, which are shown in Fig. 4 . In particular, equilibrium point E 5 exists when 15.98 < c < 201.365 and is unstable when 15.98 < c < 184. The system has periodic solutions for 15.98 < c < 184.

Fig. 4

Bifurcation diagram of equilibria with the change of c. The red, blue, pink and green curves correspond to the equilibrium points E1, E3, E4 and E5, respectively. The solid (dotted) lines indicate that the corresponding equilibrium points are stable (unstable), and H represents the Hopf bifurcation point.

Application to COVID-19

In this section, we use the mathematical model and the published data of COVID-19 spread in China to show how to deploy the vaccination program to establish an immune barrier for disease control. The inactivated vaccine BIBP developed by Sinopharm group is selected, which has 79 % protective efficacy [31] and protection duration of six months [32]. According to the medical payment regulations for confirmed and suspected patients formulated by the National Healthcare Security Administration [33], the payoff is set at 5000 Chinese Yuan. In addition, we set m = 0.4 to mean that individuals pay more attention to the daily new cases. Moreover, we fix ρ = 0.01 according to the supply of vaccine in the initial stage of vaccination in China [34]. The other parameter values in the simulations below are shown in Table 2 .

Table 2

The parameter values for the simulations of model (2.7).

Parameters	Meaning	Default values	Source
N	Total number	1.41 × 109	[35]
σ	Vaccine efficacy	0.79	[31]
1μ	Average life	75 years	[35]
β	Transmission rate	0.45 day−1	[36], [37]
1γ	Recovery time after quarantine	7 days	[38]
1γ_2	Average hospitalization days	10 days	[39]
q0	Quarantine rate	0.6 day−1	See text
ϕ	Waning rate of vaccine	0.0055 day−1	[32]
κ	Waning rate of immunity in restorers	0.0042 day−1	[40]
τ	Information interval	5 days	See text
ν	Imitation rate	0.002 day−1	See text
α	Memory decay parameter	0.03 day−1	See text
m	Weight coefficient	0.4	See text
cv	Cost of vaccination risk	500 Yuan	See text
ci	Cost of infection risk	5000 Yuan	See text
ρ	Fraction of vaccine shots	0.01 day−1	See text

Vaccination schemes for immune barrier

We focus on the immune barrier which means that all individuals in the population are protected by immunity when the vaccine coverage in the population reaches a threshold proportion. Clearly, the immune barrier is formed when the effective reproduction number R 0 ≤ 1. Using the parameter values in Table 2 and fixing θ = 0.25, numerical computations lead to R 0 = 2.2179, R 0 = 1.6348 > 1, R 1 = 0.9139 < 1. Thus, the endemic equilibrium point E 3 exists, which has no voluntary vaccination. In order to drive R 0 below 1 to suppress the endemic state, we vary θ to see how the compulsory vaccination proportion affects the reproduction number, which is shown in Fig. 5 . The numerical computations indicate R 0 ≤ 1 for θ ≥ 0.8535. Hence, the compulsory vaccination proportion must be above 85.35 % to eradicate the disease prevalence.

Fig. 5

Graph of effective reproduction number R0 versus the proportion of mandatory vaccination θ.

However, the compulsory vaccination threshold seems too high in practice. Thus, we turn to consider the importance of voluntary vaccination. Since R = 1.0001, in view of the theoretical analysis, we see that the endemic equilibrium point E 5 exists because

Assume that the vaccination cost c is so small that one neglects it to take c = 0. Then we get from (3.4), (3.5) that the infection proportion i 5 = 0 and the effective vaccination ratio is given by v 5 = 1 − 1/R 0, which means that the voluntary vaccination ratio x 5 is crucial for the herd immunity threshold. Next, we consider how the vaccination ratio x 5 for the herd immunity is influenced by the proportion of compulsory vaccination in Fig. 6 , which demonstrates that x 5 decreases with the increase of θ. More specifically, when θ = 0.6 we have x 5 = 0.6339, that is, when the proportion of compulsory vaccination is 60 %, the proportion 63.39 % of voluntary vaccination is enough to form an immune barrier. To gain intuitions how the combinations of compulsory vaccination and voluntary vaccination form the immune barrier, we present Table 3 :

Fig. 6

The proportion of voluntary vaccination x5 required to form an immune barrier is a decreasing function of the compulsory vaccination fraction θ.

Table 3

The proportions of compulsory vaccination θ and voluntary vaccination x5 required to achieve herd immunity. Vaccination coverage in each group corresponds to the herd immunity threshold calculated above.

θ	0.6	0.5	0.4	0.3	0.2
x5	0.6339	0.7071	0.7559	0.7908	0.8169
Vaccination coverage	0.85356	0.85355	0.85354	0.85356	0.85352

For the case where the perceived vaccination cost c is significant or the infection cost is small, one would expect that the individual's willingness for vaccination declines and therefore, more people don't confer enough immunity and get infected. Indeed, this is supported by the numerical simulations in panel (a) of Fig. 7 that indicates the infection proportion at equilibrium point E 5 is an increasing function of vaccination cost c . Similarly, since the weight coefficient m is inclined to the vaccination cost, the infection proportion at equilibrium point E 5 is an increasing function of m, as shown in panel (b) of Fig. 7.

Fig. 7

Influences of vaccination cost c and strategic weight coefficient m on the infection proportion i5 at equilibrium E5.

Vaccination strategies for variants

Virus variants emerge during the pandemic progression of COVID-19. These variants weaken the protective immunity of vaccines and are more contagious than the original SARS-CoV-2 strain. Table 4 shows the resultant ratios by the main variants of coronavirus in disease transmissions and damage to immunity. Fig. 8 demonstrates how these variants affects the dynamics of vaccination and infectious transmissions.

Table 4

Transmission risks and reduced efficacy of strains of coronavirus.

Strain	Found site	Incidence	Vaccine efficacy (name of vaccine)
Alpha	UK	+43 % [41]	70.4 % (ChAdOx1 nCoV-19 vaccine) [46]
Beta	South Africa	+50 % [42]	57 % (NVX-CoV2373 vaccine) [47]
Gamma	Brazil	+70 % [43]	50.38 % (Corona Vac vaccine) [47]
Delta	India	+97 % [44]	59 % (SARS-CoV-2 inactivated vaccine) [48]
Omicron	South Africa	+170.87 % [45]	65.5 %(BNT162b2 vaccine) [49]

Fig. 8

The proportions of infection for the five mutant strains as time evolves.

The epidemic Omicron strain is selected as an example to analyze what kind of vaccine strategy should be adopted for mutant strain. For the Omicron mutant strain, we obtain R 0 = 6.0076. Considering that the supply of vaccine during the spread of Omicron virus is much increased, we set ρ = 0.03. For the best case where all individuals are vaccinated with two doses of Pfizer vaccine, i.e., θ = 1, the effective reproduction number R 0 = R 1 = 1.3206 > 1. This means that the vaccine is not effective enough to provide protection against the transmission of Omicron strain. Note that the efficacy of the vaccine wanes over time and booster shots improve immune protection. In fact, the data in the UK shows that the maximum effectiveness of booster dose in preventing symptoms is up to 75 % [49], [50], which is even better than the best vaccine efficacy of first two doses. In order to observe the relation between the vaccine efficacy and the corresponding effective reproduction number, we present Table 5 , which indicates that the effective reproduction number becomes larger as the efficacy of the vaccine declines. More importantly, the facts that every reproduction number in Table 5 is larger than unity means the vaccine would not provide adequate protection against omicron virus, even though full vaccinations are practiced because of θ = 1.

Table 5

The effective reproduction number R0 is always larger than unity and is increased with the decrease of vaccine efficacy σ where θ = 1.

σ	0.75	0.6	0.65	0.6	0.55	0.5
Rv0	1.1864	1.2534	1.3285	1.4132	1.5094	1.6196

Since the vaccinations are insufficient to combat the spread of Omicron virus, we are forced to speed up the quarantine rate. Fig. 9(a) demonstrates that the effective reproduction number R 0 is a decreasing function of the quarantine rate q 0 when σ = 0.75 and θ = 1. More specifically, we have R 0 ≤ 1 for q 0 ≥ 0.978. We point out that the constraint θ = 1 is not necessary. Indeed, Fig. 9(b) indicates that a larger quarantine rate q 0 leads to a reasonable compulsory ratio θ for R 0 = 1.

Fig. 9

Panel (a) shows that the effective reproduction number R0 is decreasing with quarantine rate q0 and R0 ≤ 1 for q0 ≥ 0.978. Panel (b) indicates that a reasonable compulsory vaccine coverage θ is applicable for an adequate quarantine rate q0 in order to form an immune barrier.

Let us fix θ = 0.25, which is a reasonable compulsory vaccination proportion of medical personnel and service workers in China. Then for reasonable quarantine rates, this proportion of compulsory vaccination cannot reach the vaccine coverage required for the formation of immune barrier. Thus, it is important to add a proportion of voluntary vaccinations. As in Table 3, we present Table 6 to show how the combinations of voluntary vaccination with quarantine rate form the immune barrier and how the vaccination coverage is needed. Clearly, a lower vaccination coverage needs a faster quarantine rate to control the outbreak of epidemic disease.

Table 6

The combinations of voluntary vaccination proportion with quarantine rate to form the immune barrier and their corresponding vaccination coverage where θ = 0.25.

Quarantine rate q0	1	1.5	2	2.5	3
Voluntary vaccination x5	0.985	0.704	0.5049	0.3564	0.2415
Vaccine coverage	0.9888	0.778	0.6286	0.5173	0.4311

Discussions

In this paper, we propose the mathematical model with both compulsory vaccination and voluntary vaccination, and use the game theory to describe the behavior dynamics of individuals who consider the cost and benefit of vaccination. To mimic the realistic decision-making process, we take the decision function that is distributed on an interval for the past information. Furthermore, we adopt the nonlinear response function to the pay-off of vaccination, which is an extension of the classical Fermi function in discrete imitation dynamics. We compute the basic reproduction number R 0 and the effective reproduction number R and analyze the existence and stability of the equilibrium points of the model. We find that the model admits the Hopf bifurcation when the imitation intensity is larger or the length of information interval is larger.

We use the model (2.7) to simulate the dynamics of COVID-19 spread in China. The influences of the key parameters on disease prevalence are given in Fig. 7. With the aid of the published data, we show how to deploy the vaccination program to establish an immune barrier for disease control. First, we compute the threshold value of vaccine coverage for immune barrier where there are only mandatory vaccinations, which looks too high to be implemented in practice. For this reason, we calculate the combinations of compulsory vaccination with voluntary vaccination to form the immune barrier in Table 3. This is of importance in designing the vaccination scheme. Indeed, once an upper percentage of the individuals with voluntary vaccination is estimated, one can select from the combination the proportion of mandatory vaccination to establish the immune barrier in practice.

Big challenges in fighting pandemic COVID-19 are the emergence of virus variants, which weaken the protective immunity of vaccines and are more contagious than the original SARS-CoV-2 strain. We show that the effective reproduction number R 0 for Omicron strain is larger than unity even though all individuals are vaccinated with two doses of Pfizer vaccine. Then we consider the effect of booster shots in Table 5, which shows that the effective reproduction numbers are always larger than unity for different vaccine efficacy and full mandatory vaccination. As a result, we strengthen the vaccination effect by superposing the quarantine of infectious individuals. Table 6 presents the combinations of voluntary vaccination proportion and quarantine rate to form the immune barrier and their corresponding vaccination coverage where the proportion of mandatory vaccination is fixed as θ = 0.25. This is helpful to select an optimal quarantine intensity on the status of population immunity in practice. In all, our mathematical modelling and analysis give the theoretical basis to design an optimal strategy of vaccination, which is aided with quarantine, to control a high contagious variant of infectious disease.

The mathematical model (2.7) captures the main features of infectious diseases with vaccination games. Notice that the latent period and asymptomatic phase of some diseases such as COVID-19 and measles contribute more or less to the disease outbreaks. It will be interesting to introduce latent compartment and asymptomatic compartment into the model to describe the disease transmissions more accurately. We leave these as future researches.

CRediT authorship contribution statement

Jingwen Ge: Investigation, Writing – original draft. Wendi Wang: Supervision, Investigation, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.