1 Introduction

Vaccines have saved many lives throughout our history (Domachowske and Suryadevara 2021). For instance, in the period 2000–2019, vaccination prevented 37 million of deaths in 98 countries (Li et al. 2021). However, despite reducing mortality, reaching a high vaccine coverage remains a challenge worldwide, as alerted by the World Health Organization (2019). Two barriers are availability and confidence. Vaccine availability depends on supply and logistic networks (Zaffran et al. 2013); vaccine confidence depends on perceived risks and benefits (Lane et al. 2018).

Sociocultural, political, and religious factors influence the likelihood of getting a vaccine (Lane et al. 2018; Lin et al. 2021; Loomba et al. 2021; Megget 2020; Neumann-Böhme et al. 2020). In addition, the diffusion of misinformation and factual information affects the public perception about safety and efficacy of vaccines (Lane et al. 2018; Lin et al. 2021; Loomba et al. 2021; Megget 2020; Neumann-Böhme et al. 2020). Unfortunately, the anti-vaccine sentiment, fed by irrational beliefs and conspiracy theories (Latkin et al. 2021), is widespread and can undermine the control of epidemics by vaccination (Hussain et al. 2018; Lane et al. 2018; Smith and Graham 2019). In fact, if the proportion of the unvaccinated group remains above a critical number, herd immunity can not be achieved and the pathogen can endemically persist in the host population (Anderson and May 1985). Hence, the anti-vaccine movement has been identified as a top threat to global health (World Health Organizations 2019). COVID-19, the contagious disease responsible for a huge public health emergency and an unprecedent economic crisis, has fueled the debate on vaccine hesitancy (Latkin et al. 2021; Lin et al. 2021; Loomba et al. 2021; Megget 2020; Neumann-Böhme et al. 2020). However, vaccination can be the only way to achieve and maintain a sustained reduction in COVID-19 infections (Priesemann et al. 2021a, b).

The vaccination decision-making process has been considered in theoretical investigations based on ordinary differential equations (ODE) (Arefin et al. 2020; Buonomo 2020; Buonomo et al. 2022; Flaig et al. 2018), game theory (Schimit and Monteiro 2011; Liu et al. 2012; Wang et al. 2016), regular and complex networks (Curiel and Ramirez 2021; Pires et al. 2018, 2021), Bayesian inference (Coelho and Codeço 2009). In these studies, the propagation of infections in human communities is affected by beliefs, fears, circulating rumors, factual information. Here, a model formulated in terms of ODE is used for studying the spread of vaccine-preventable infectious diseases, by taking into account the vaccination intention of the susceptible individuals. Thus, each susceptible individual is either pro-vaccine or anti-vaccine; however, change of mind is possible.

The text is organized as follows. The proposed model is introduced in Sect. 2. Analytical results are presented in Sect. 3 and numerical simulations illustrating these results are shown in Sect. 4. In Sect. 5, the possible relevance of this work is stressed by considering the current COVID-19 outbreak.

2 The Epidemic Model

Let the host population be divided into four groups: the P-group of the susceptible individuals who are pro-vaccine, the A-group of the susceptible individuals who are anti-vaccine, the I-group of the infected individuals; and the R-group of the recovery individuals. Assume that these four groups are uniformly distributed over the region where they live; thus, the time evolution of the number of individuals belonging to each group is not affected by the particular geographic location of its individuals. In other words, assume that the homogeneous mixing assumption holds (Turnes and Monteiro 2014). The PAIR model proposed here is written as the following set of ODE:

$$\begin{aligned} \frac{dP(t)}{dt}= & {} -aP(t)I(t) + f[cI(t) + hR(t)] - vP(t) + \alpha A(t) - \beta P(t) \end{aligned}$$
(1)
$$\begin{aligned} \frac{dA(t)}{dt}= & {} -aA(t)I(t) + (1-f)[cI(t) + hR(t)] - \alpha A(t) + \beta P(t) \end{aligned}$$
(2)
$$\begin{aligned} \frac{dI(t)}{dt}= & {} a[P(t)+A(t)]I(t) -bI(t)-cI(t)\end{aligned}$$
(3)
$$\begin{aligned} \frac{dR(t)}{dt}= & {} bI(t)- hR(t) + vP(t) \end{aligned}$$
(4)

In these equations, the variables P(t), A(t), I(t), and R(t) denote the numbers of P, A, I, and R-individuals at the time t, respectively; a is the infection rate constant of P and A-individuals; b the recovery rate constant of I-individuals; c the death rate constant of I-individuals; h the death rate constant of R-individuals; v the immunization rate constant (that is, the product between the vaccination rate constant and the vaccine efficacy) of P-individuals; \(\alpha\) the rate constant representing migration from the A-group to the P-group; and \(\beta\) the rate constant representing migration from the P-group to the A-group. Disseminating negative false news about vaccines can decrease \(\alpha\) and increase \(\beta\); disseminating data related to the protective power of vaccines can increase \(\alpha\) and decrease \(\beta\). Thus, the influences of traditional media and social media are taken into account in the parameters \(\alpha\) and \(\beta\).

Notice that the contagious disease is transmitted through contacts between P and I-individuals and between A and I-individuals according to the same rate constant a. Thus, it is assumed that the P-group and the A-group take similar preventive measures (such as frequent hand washing, wearing face masks, use of latex condoms). Also, it is assumed that both groups have similar spatial movement patterns and are aware that the pathogen can cause a serious illness. The difference between them is that P-individuals trust in the benefits of getting a vaccine and A-individuals do not.

Notice also that the deaths of I and R-individuals are balanced by the births of P and A-individuals. In 2020, the growing rate (per year) of the population was, for instance, 0.59% in USA and 0.04% in Russia (Worldometer 2021). Thus, supposing that the deaths are balanced by the births can be considered a reasonable assumption. Observe that this assumption implies that the birth of individuals that remain susceptible throughout their lives is balanced by the death of these individuals at the same rate constant. Since these terms cancel out, they were omitted in Eqs. (1) and (2). In addition, in this model, individuals cannot be born infected or immunized.

The parameter f denotes the fraction of births of P-individuals (the fraction of susceptible individuals who engage in the P-group); consequently, \(1-f\) is the fraction of births of A-individuals. When R-individuals lose their acquired immunity, they return to the susceptible groups according to these same proportions. In this model, the parameters a, b, c, f, h, v, \(\alpha\), and \(\beta\) are non-negative numbers.

Since \(dP(t)/dt+dA(t)/dt+dI(t)/dt+dR(t)/dt=0\), then \(P(t)+A(t)+I(t)+R(t)=N\); that is, the total number of individuals remains constant and equal to N. Therefore, \(R(t)=N-P(t)-A(t)-I(t)\) and the original set of ODE can be rewritten as:

$$\begin{aligned} \frac{dP}{dt}= & {} -aPI + f[cI + h(N-P-A-I)] - vP + \alpha A - \beta P \equiv F_1(P,A,I) \end{aligned}$$
(5)
$$\begin{aligned} \frac{dA}{dt}= & {} -aAI + (1-f)[cI + h(N-P-A-I)] - \alpha A + \beta P \equiv F_2(P,A,I) \end{aligned}$$
(6)
$$\begin{aligned} \frac{dI}{dt}= & {} a(P+A)I -bI-cI \equiv F_3(P,A,I) \end{aligned}$$
(7)

From a public health point of view, the best scenario is given by \(f=1\) (only P-individuals can be introduced in the population) and \(\beta =0\) (no migration to the A-group), and the worst scenario is given by \(f=0\) (only A-individuals can be introduced in the population) and \(\alpha = 0\) (no migration to the P-group).

In the next section, this model is analyzed from a dynamical systems theory perspective (Guckenheimer and Holmes 2002).

3 Analytical Results

In the state space \(P \times A \times I\), an equilibrium point \((P^*,A^*,I^*)\) corresponds to a stationary solution; that is, a solution satisfying \(dP(t)/dt=0\), \(dA(t)/dt=0\), and \(dI(t)/dt=0\). Such a solution is given by \(P(t)=P^*\), \(A(t)=A^*\), and \(I(t)=I^*\), in which the constants \(P^*\), \(A^*\), and \(I^*\) are obtained from \(F_1(P^*,A^*,I^*)=0\), \(F_2(P^*,A^*,I^*)=0\), and \(F_3(P^*,A^*,I^*)=0\).

The third-order dynamical system given by Eqs. (5)–(7) presents a disease-free equilibrium point (that is, an equilibrium point with \(I^*=0\)). Its coordinates are:

$$\begin{aligned} P_{free}^*= & {} \frac{\alpha h N}{[\alpha + \beta + (1-f)v]h + \alpha v} \end{aligned}$$
(8)
$$\begin{aligned} A_{free}^*= & {} \frac{[\beta +(1-f)v] h N}{[\alpha + \beta + (1-f)v]h + \alpha v} \end{aligned}$$
(9)
$$\begin{aligned} I_{free}^*= & {} 0 \end{aligned}$$
(10)

Obviously, \(R_{free}^* = N - P_{free}^* - A_{free}^* - I_{free}^*\).

The local stability of an equilibrium point \((P^*,A^*,I^*)\) can be inferred from the eigenvalues \(\lambda\) of the Jacobian matrix J, which corresponds to the linearization of Eqs. (5)–(7) around this point (Guckenheimer and Holmes 2002). Here:

$$\begin{aligned}&{{\textbf{J}}}(P^{*},A^{*},I^{*}) \;\nonumber \\&= \; \left[ \begin{array}{ccc} -aI^* - fh - v - \beta &{} -fh + \alpha &{} -aP^* + f(c-h) \\ -(1-f)h + \beta &{} -aI^* - (1-f)h - \alpha &{} -aA^* + (1-f)(c-h) \\ aI^* &{} aI^* &{} a(P^*+A^*) - b - c\\ \end{array} \right] \end{aligned}$$
(11)

An equilibrium point is locally asymptotically stable if the real part of all eigenvalues \(\lambda\) of J is negative (Guckenheimer and Holmes 2002). Recall that \(\lambda\) is determined from \(\det (\mathbf{{J}}- \lambda \mathbf{I})=0\) (in which \(\mathbf{I}\) is the identity matrix). Thus, the eigenvalues \(\lambda\) associated with the disease-free solution are the roots the polynomial:

$$\begin{aligned} (\lambda ^{2} + \theta _{1} \lambda + \theta _{2}) [a(P^{*}+A^{*}) - b - c - \lambda ]=0 \end{aligned}$$
(12)

with \(\theta _1 =\alpha +\beta +h+v\) and \(\theta _2=[\alpha + \beta + (1-f)v]h + \alpha v\). Since \(\theta _1 >0\) and \(\theta _2 > 0\), then both roots of \(\lambda ^2 + \theta _1 \lambda + \theta _2=0\) have negative real part. Therefore, the disease-free equilibrium point is locally asymptotically stable if the third eigenvalue is negative. This condition can be expressed as \(R_0 < 1\) if \(R_0\) is defined as:

$$\begin{aligned} R_{0} \equiv \frac{a(P_{free}^{*} + A_{free}^{*})}{b+c} = \left( \frac{aN}{b+c} \right) \frac{[\alpha + \beta + (1-f)v]h }{[\alpha + \beta + (1-f)v]h + \alpha v} \end{aligned}$$
(13)

In epidemic models, \(R_0\) is a parameter known as basic reproduction number (Anderson and May 1992; Schimit and Monteiro 2012). Usually, for \(R_0 < 1\), the disease-free equilibrium point is asymptotically stable and an endemic equilibrium point (that is, an equilibrium point with \(I^*>0\)) is unstable; for \(R_0>1\), these equilibrium points exchange their stabilities. Therefore, \(R_0\) is a bifurcation parameter, because the corresponding dynamical system experiences a (transcritical) bifurcation (Guckenheimer and Holmes 2002) for a particular value of \(R_0\) (in this case, this particular value is \(R_0=1\)). This parameter is also interpreted as the average number of secondary infections caused by a single infected individual inserted into a completely susceptible population (Anderson and May 1992; Schimit and Monteiro 2012). Hence, if \(R_0 < 1\), the disease tends to naturally disappear; if \(R_0>1\), the disease tends to endemically persist.

For \(\alpha =0\) and \(\beta =0\), then \(R_0 = r_0\), with:

$$\begin{aligned} r_{0} \equiv \frac{aN}{b+c} \end{aligned}$$
(14)

This expression was derived in other works on theoretical epidemiology (Ferraz and Monteiro 2019; Schimit and Monteiro 2012). Note that \(R_0 \le r_0\) for \(\alpha \ge 0\), \(\beta \ge 0\), and \(v \ge 0\). Thus, \(r_0\) is the superior limit of \(R_0\).

Equation (13) can be used to derive the critical immunization rate constant \(v_{c}\) above which the disease can be eradicated by vaccination. Observe that \(R_0 < 1\) implies:

$$\begin{aligned} v > v_{c} \equiv \frac{(\alpha +\beta )h(r_0-1)}{(1-f)h(1-r_0)+\alpha } \end{aligned}$$
(15)

For \(R_0>1\), the system converges to the endemic equilibrium point given by:

$$\begin{aligned} P_{ende}^*= & {} q - pI_{ende^*} \end{aligned}$$
(16)
$$\begin{aligned} A_{ende}^*= & {} \frac{N}{r_0} - P_{ende^*} \end{aligned}$$
(17)
$$\begin{aligned} I_{ende}^*= & {} \frac{-y+\sqrt{y^2-4xz}}{2x} \end{aligned}$$
(18)

in which \(p=(b+h)/v\), \(q=hN(r_0-1)/(v r_0)\), \(x=ap\), \(y=-aq+f(c-h)+p(\alpha +\beta +v)\), and \(z=(\alpha N/r_0)-q[\alpha +\beta +(1-f)v]\), with \(z<0\) (that is, \(v<v_c\)). Obviously, \(R_{ende}^* = N - P_{ende}^* - A_{ende}^* - I_{ende}^*\).

In the best scenario (that is, \(f=1\) and \(\beta =0\)), then \(P_{ende}^*=N/r_0\), \(A_{ende}^*=0\), and \(I_{ende}^*=hN[1-(1/R_0)]/(b+h)\), with \(R_0 = r_0 h/(h+v)\). Therefore, the higher the value of v, the smaller the values of \(R_0\) and \(I_{ende}^*\). Also, in this scenario, the A-group vanishes. In the worst scenario (that is, \(f=0\) and \(\alpha =0\)), then \(P_{ende}^*=0\), \(A_{ende}^*=N/r_0\), and \(I_{ende}^*=hN[1-(1/r_0)]/(b+h)\). Since \(R_0=r_0\), then \(R_0\) and \(I_{ende}^*\) can not be reduced by taking \(v>0\). In addition, the P-group disappears. Note that, except in these two extreme scenarios, the exact formulas for the coordinates of the endemic equilibrium point given by Eqs. (16)–(18) are so intricate that approximate formulas may be more useful in practice.

Consider that \(v \ll 1\) (this choice is justified in the next section). In this limit, approximate expressions for the attractor endemic steady-state are given by:

$$\begin{aligned} P_{approx}^*= & {} \frac{N}{r_0}-A_{approx}^{*} \end{aligned}$$
(19)
$$\begin{aligned} A_{approx}^*\approx & {} \frac{(\beta N/r_0)+[(1-f)hN(r_0-1)/r_0]+(1-f)(c-h)I_{approx}^*}{\alpha +\beta + aI_{approx}^{*}} \end{aligned}$$
(20)
$$\begin{aligned} I_{approx}^*\approx & {} \frac{hN}{(b+h)} \left( 1 - \frac{1}{r_0} \right) \end{aligned}$$
(21)

Observe that Eq. (21) (the approximate endemic solution obtained by taking \(v \ll 1\)) is equal to Eq. (18) for the worst scenario (the exact endemic solution obtained by taking \(f=0\) and \(\alpha =0\)). Therefore, in this context, low immunization rate is equivalent to no influx to the pro-vaccine group.

Observe also that, for \(v \ll 1\), \(I_{approx}^*\) does not depend on \(\alpha\), \(\beta\), and f. Therefore, these three parameters do not affect the long-term amount of sick people (that is, \(I_{approx}^*\)). They do not affect also the long-term amount of recovered people, because \(R_{approx}^* = N - (N/r_0)-I_{approx}^*\). However, these three parameters affect the long-term composition of the susceptible population (that is, \(P_{approx}^*\) and \(A_{approx}^*\)) and the critical immunization rate constant \(v_c\) given by Eq. (15).

In the next section, these analytical results are illustrated by computer simulations performed with plausible parameter values for COVID-19. The propagation of this contagious disease has been analytically and numerically studied by using similar epidemiological models (Angeli et al. 2022; Grimm et al. 2021; Harari and Monteiro 2021; He et al. 2020; Ram and Schaposnik 2021).

4 Numerical Results

Equations (1)–(4) were numerically solved by employing the 4th-order Runge-Kutta integration method with integration time step of 0.01. In the simulations, \(N=1\); thus, the variables P(t), A(t), I(t), and R(t) correspond to the percentages of P, A, I, and R-individuals in the whole population. The time t is measured in weeks. The initial condition of the simulations is \(P(0)=0.49\), \(A(0)=0.49\), \(I(0)=0.02\), and \(R(0)=0\); however, as in other works (Monteiro et al. 2020), the long-term behavior does not depend on the initial condition (it depends only on the parameter values).

Here, parameter values related to the infection by SARS-CoV-2 are considered. Assume that the average life expectancy of the population is 80 years; assume also that the immunity (either acquired after recovery or elicited by vaccine) lasts for 1 year, on average. Therefore, \(h=1/(80 \times 52) + 1/52\). Assume that the recovery period is about 2 weeks (Gostic et al. 2020; Singhal 2020); thus, \(b=1/2\). In addition, assume that I-individuals can die only due to the infection. Consider that 3% of I-individuals die in 3 weeks (Johns Hopkins University 2021; Singhal 2020); hence, \(c=1/100\). Also, take \(f=3/4\). This choice was motivated by surveys on vaccine acceptance. For instance, surveys conducted in 7 European countries showed that 73.9% of the respondents are willing to get vaccinated against COVID-19 (Neumann-Böhme et al. 2020). This percentage is lower in USA (Loomba et al. 2021; Perlis et al. 2020) and higher in China (Chen et al. 2021). With this value of f, the aim is to obtain the fraction of the susceptible population who is pro-vaccine similar to those found in such surveys.

In addition, take \(a=3/5\). With this choice, for \(\alpha =0\), \(\beta =0\), and \(v=0\), the basic reproduction number of COVID-19 is \(R_0 = r_0 \simeq 1.2\). A similar value was found in the beginning of the pandemic after imposing travel restrictions (Kucharski et al. 2020) and is still found in many countries (EpiForecasts 2021). In the absence of preventive measures, \(R_0\) can be higher (Gostic et al. 2020; Singhal 2020).

Between the end of 2020 and the beginning of 2021, the USA had fully vaccinated 25 million people in 10 weeks (Our World in Data 2021). Since its population is 330 million people, then \(v = 0.9 \times 25/(330 \times 10) \simeq 7/1000\), by considering 90% vaccine efficacy (Kim et al. 2021). In the UK, 0.8 out of 67 million people were fully vaccinated in 11 weeks; therefore, \(v = 0.9 \times 0.8/(67 \times 11) \simeq 1/1000\). In many other countries, however, v has been lower (Our World in Data 2021). Hence, here, \(v = 0.5/1000\).

Figure 1a, b show the time evolutions of P(t) (blue line), A(t) (cyan line), I(t) (red line), and R(t) (green line). In Fig. 1a, \(\alpha = 1/10\) and \(\beta = 1/20\); in Fig. 1b, \(\alpha = 1/20\) and \(\beta = 1/10\). In both cases, \(R_0 \simeq 1.2\). In Fig. 1a, \(P(t) \rightarrow 0.568\), \(A(t) \rightarrow 0.282\), \(I(t) \rightarrow 0.005\), and \(R(t) \rightarrow 0.145\), as time t progresses. These numbers are in good agreement with the values of \(P^*\), \(A^*\), \(I^*\), and \(R^*=1-(P^*+A^*+I^*)\) obtained from the analytical expressions given by Eqs. (20)–(21). In Fig. 1b, \(P(t) \rightarrow 0.291\), \(A(t) \rightarrow 0.559\), \(I(t) \rightarrow 0.005\), and \(R(t) \rightarrow 0.145\). In this case, the coordinates of the endemic equilibrium point \((P^*,A^*,I^*)\) reached in the numerical simulation also agree with the numbers computed from Eqs. (20)–(21).

In the proposed model, the higher the value of \(\beta\), the greater the spread of misinformation about vaccines; the higher the value of \(\alpha\), the greater the effort to combat these fake news. Figure 1a, b show that, as expected, \(P^*\) increases with \(\alpha\) and \(A^*\) increases with \(\beta\). In Fig. 1a, \(P^*/(P^*+A^*) \simeq 67\%\), which is the fraction of the susceptible population who is pro-vaccine found in the USA (Lin et al. 2021; Perlis et al. 2020). In Fig. 1b, this fraction is \(34\%\).

For the parameter values of Fig. 1a, disease eradication is achieved by vaccination for \(v > v_c \simeq 5.3/1000\); in Fig. 1b, for \(v > v_c \simeq 11/1000\). Notice that the estimated value of v for the USA (7/1000) in the beginning of 2021 is enough to eliminate the disease in the scenario corresponding to Fig. 1a, but is not enough in the scenario of Fig. 1b. Hence, reducing misinformation about vaccines can facilitate disease eradication. In the end of March of 2022, there were 2175 million people fully immunized in USA (Johns Hopkins University 2022); therefore, \(v \simeq 9/1000\). This value is still not enough to eradicate the disease in the second scenario.

Figure 2a, b exhibit the time evolutions of P(t) (blue line), A(t) (cyan line), I(t) (red line), and R(t) (green line) for \(v=15/1000\), which surpasses the critical rate \(v_c\) in both plots. In these figures, the other parameter values are the same as those used in Fig. 1a, b, respectively. In Fig. 2a, \(P(t) \rightarrow 0.433\), \(A(t) \rightarrow 0.233\), \(I(t) \rightarrow 0\), and \(R(t) \rightarrow 0.334\); in Fig. 2b, \(P(t) \rightarrow 0.260\), \(A(t) \rightarrow 0.540\), \(I(t) \rightarrow 0\), and \(R(t) \rightarrow 0.200\). In both cases, the coordinates of the disease-free equilibrium points reached in the numerical simulations are equal to the numbers computed from Eqs. (8)–(10).

Fig. 1
figure 1

Numerical solutions of Eqs. (1)–(4) for \(a=3/5\), \(b=1/2\), \(c=1/100\), \(f=3/4\), \(h=1/(80 \times 52) + 1/52\), and \(v=0.5/1000\). In (a), \(\alpha = 1/10\) and \(\beta = 1/20\); in (b), \(\alpha = 1/20\) and \(\beta = 1/10\). Also, \(N=1\); hence, P(t) (blue line), A(t) (cyan line), I(t) (red line), and R(t) (green line) represent normalized concentrations. The time t is measured in weeks. The initial condition is \(P(0)=49\%\), \(A(0)=49\%\), \(I(0)=2\%\), and \(R(0)=0\). In both plots, \(I(t) \rightarrow 0.5\%\), \(R(t) \rightarrow 14.5\%\), and \(P(t)+A(t) \rightarrow 85\%\). In (a), \(P(t) \rightarrow 57\%\) and \(A(t) \rightarrow 28\%\); in (b), \(P(t) \rightarrow 29\%\) and \(A(t) \rightarrow 56\%\). Notice that, for \(v \ll 1\), the parameters \(\alpha\) and \(\beta\) do not influence the asymptotical value of the sum \(P(t)+A(t)\), but they influence the asymptotical value of P(t) and the asymptotical value of A(t). Also, notice that P(t) and A(t) can present equivalent dynamical behaviors if the values of \(\alpha\) and \(\beta\) are switched as done in (a) and (b). (Color figure online)

Full size image
Fig. 2
figure 2

Numerical solutions of Eqs. (1)–(4) for \(a=3/5\), \(b=1/2\), \(c=1/100\), \(f=3/4\), \(h=1/(80 \times 52) + 1/52\), and \(v=15/1000\). In (a), \(\alpha = 1/10\) and \(\beta = 1/20\); in (b), \(\alpha = 1/20\) and \(\beta = 1/10\). In addition, \(N=1\); therefore, P(t) (blue line), A(t) (cyan line), I(t) (red line), and R(t) (green line) represent normalized concentrations. The time t is measured in weeks. The initial condition is \(P(0)=49\%\), \(A(0)=49\%\), \(I(0)=2\%\), and \(R(0)=0\). In (a), \(P(t) \rightarrow 43\%\), \(A(t) \rightarrow 23\%\), \(I(t) \rightarrow 0\), and \(R(t) \rightarrow 34\%\); in (b), \(P(t) \rightarrow 26\%\), \(A(t) \rightarrow 54\%\), \(I(t) \rightarrow 0\), and \(R(t) \rightarrow 20\%\). In both cases, the disease is eradicated by vaccination because \(v > v_c\). (Color figure online)  

Full size image

5 Discussion and Conclusion

Here, an epidemic model based on ODE was proposed and analyzed. This model describes the transmission dynamics of a vaccine-preventable infection, by taking into consideration pro and anti-vaccination groups. The analytical expressions derived in Sect. 3 predict the long-term composition of a population subject to a vaccine-preventable infection that endemically persists, as measles and possibly COVID-19. In fact, COVID-19 is still a major problem (Allaerts 2021) and it can remain in our population after the pandemic is over (John and Seshadri 2021).

It was shown that if \(R_0<1\), the disease-free steady-state, given by Eqs. (8)–(10), is reached; if \(R_0>1\), the attractor is the endemic steady-state, given by Eqs. (16)–(18). For \(v \ll 1\), the coordinates of the endemic attractor can be approximated by Eqs. (20)–(21). These approximated expressions reveal that, for \(v \ll 1\), the values of the parameters \(\alpha\), \(\beta\), f, and v do not significantly affect \(I^*\), \(R^*\), and \(P^*+A^*\); they only affect \(P^*\) and \(A^*\); that is, the composition of the susceptible group. Therefore, if \(v \ll 1\), the fight against false information about vaccines may not have a meaningful impact on the long-term number of sick individuals. This finding is illustrated by Fig. 1a, b, which were plotted by taking parameter values of the current COVID-19 pandemic. In these figures, \(I(t) \rightarrow I_{approx}^* \simeq 0.5\%\). In the USA, this fraction corresponds to 1.65 million people infected by SARS-CoV-2 at any time step. Recall that this number was obtained by taking \(v=0.5/1000\), which is a low value of v. In fact, the current vaccination rate in this country is \(v \simeq 9/1000\). With this value, \(I(t) \rightarrow I_{free}^* = 0\) in the scenario of Fig. 1a and \(I \rightarrow I_{ende}^* \simeq 0.077 \%\) in the scenario of Fig. 1b (about 250 thousand people infected in the USA at any time step). It is relevant to stress that any prediction for \(I^*\) can be tested only after I(t) reaching its long-term behavior and it depends on v. The available data, however, still reflect the transient phase of this disease spread.

Notice that the condition \(v \ll 1\) can describe the vaccination rate of other contagious diseases, such as measles, mumps, and rubella, because only children usually take MMR vaccine (Kauffmann et al. 2021), and hepatitis A, due to the poor adherence to the vaccination against this infection (Johnson et al. 2019).

Although \(\alpha\), \(\beta\), and f do not significantly influence the long-term number of sick individuals for \(v \ll 1\) , these three parameters do affect the critical immunization rate constant \(v_c\). Equation (15) shows that, for instance, \(v_c\) increases with \(\beta\) (the parameter related to the transition \(P \rightarrow A\)). Thus, the reduction in vaccine acceptance found in surveys conducted in the USA in 2020 (Lin et al. 2021) implies an increase in the immunization rate required to eliminate the pathogen causing COVID-19. Such a reduction may continue to occur for the next years (Johnson et al. 2020). Hence, campaigns against fake news about vaccines are crucial for reducing (\(\beta\) and consequently) \(v_c\). Figure 2a, b confirm that, if \(v > v_c\), vaccination eradicates the contagious disease. Therefore, eradication is possible even if A-individuals remain unvaccinated.

This study predicts the long-term behavior of an epidemic model with pro and anti-vaccine groups. For COVID-19, the accuracy of the predictions can be affected by several factors, such as implementation of irregular quarantine/lockdown regimes, emergence of variants of the virus SARS-CoV-2, the use of different vaccines with different immunization rates, data used in the predictions were obtained in the transient phase of this disease propagation. However, the main message here is: to contain the spread of vaccine-preventable infections, nations must promote high vaccination engagement by mitigating concerns about effectiveness and side-effects. Misinformation hampers eradication.