Understanding the oscillations of an epidemic due to vaccine hesitancy

1.   Introduction

Vaccination is one of the most cost-effective ways of avoiding disease—it prevents 2–3 million deaths a year, and a further 1.5 million could be avoided if global coverage of vaccinations improved ^[1]. The recent COVID-19 pandemic, caused by SARS-CoV-2, presented an unprecedented challenge to human health, the economy, and nearly all aspects of our society. Vaccination, social distancing, and mask-wearing are effective mitigation methods to prevent widespread morbidity and mortality from the infection. The rapid development of effective vaccinations was one of the longstanding legacies of the COVID-19 pandemic.

Mass vaccination for COVID-19 started on December 14, 2020 in the United States ^[2], leading to a significant reduction of infections and disease burden within vaccinated populations ^[3]. As the vaccination programs progressed into March 2021, a widespread reduction in cases and deaths ^[4] led to the relaxation of non-pharmaceutical interventions in many places ^[5].

However, vaccine hesitancy, a strong unwillingness to adopt vaccination ^[6,7], remains a persistent barrier to complete population inoculation against highly infectious diseases. "It is difficult to comprehend that some of the greatest triumphs of medical science are being eroded by misinformation and distrust." ^[8]. As COVID-19 vaccination efforts went underway worldwide, significant vaccine hesitancy and resistance were observed even among healthcare workers ^[9,10]. Multi-country surveys have shown wide country-to-country variation in vaccine acceptance (31–89%) ^[11,12]. Another international study using data from 13 countries showed a concerning trend: As the pandemic progressed, the percentage of people intending to refuse vaccination increased ^[13]. In the UK, a large survey of 5114 adults in October 2020 showed 16.6% were "very unsure", and 11.7% were "strongly hesitant" ^[14]. An excellent model fit explaining 86% of the variance in hesitancy "was provided by beliefs about the collective importance, efficacy, side-effects, and speed of development of a COVID-19 vaccine."^[14] In Japan, a country with more than 75% vaccination rate, nearly 4% of the people surveyed switched from being willing to be vaccinated to being hesitant between February 2021 and February 2022 ^[15]. A more recent analysis of Polish COVID-19 vaccination data (until January 2023) showed a general decline in subsequent vaccine uptake even after increased government-provided doses, which led them to recommend significant changes in official vaccination policies ^[16].

There are many possible reasons for vaccine hesitancy ^[17], including concerns about the safety and side effects of the vaccine. Some studies suggest that belief in conspiracy theories drives vaccination hesitancy through exposure to misinformation ^[18,19]. Even medically informed individuals rely on a personal risk–benefit perception that may be influenced by misinformation regarding vaccine safety ^[9]. Distrust in the scientific expertise and health and government authorities drives consumers from "traditional" sources (newspapers, television, radio, government agencies) to social media outlets ^[20], who are then more likely subject to misinformation. Misinformation has been shown to correlate with vaccine hesitancy, even though the causal relationship has yet to be established ^[18,19].

Our goal is to develop a model to help understand vaccine hesitancy's mechanistic role in the spread of infections. Meant to explore the essentials of a system that can accommodate both stable endemic equilibria and periodic oscillations of high/low infections, our model does not address the underlying causes for vaccine hesitancy, nor attempt to reproduce the COVID-19 pandemic spread quantitatively. The key ingredient here is an evolutionary game theoretical framework for aggregate decision-making on vaccination with a susceptible-infected-susceptible ($SIV$) model for epidemic spread, which we call the $SIV$ model. We show analytically and numerically that the competition between the fear of vaccination and the perceived risk of infection introduced in our model, may drive the system into a stable endemic equilibrium or dynamic oscillations of the epidemic.

2.   A SIV model for vaccine hesitancy

The simplest epidemic model considers a well-mixed population, a fraction of which are susceptible ($S$) to the infection, and the rest being infected and infectious ($I = 1-S$). The transmission rate is $\beta$, and the infected recover at the rate $\gamma$. Known as the $SIV$ model, this system evolves in time ($t$) according to the equation: $\dot{S}\; = \; dS/dt\; = \; \gamma I-\beta IS$ (or $\dot{I} = \beta SI-\gamma I$). In general, if $\beta > \gamma$ (infection faster than recovery), the system settles into the disease endemic state with a fraction of susceptible at $\gamma /\beta$ and infectious at $1 - \gamma/\beta$. Otherwise, the system eventually arrives at a disease-free, stable, steady state with $I = 0$. This disease-free steady state loses stability when $\beta > \gamma$ or $R_0 = \beta/\gamma > 1$. The $SIV$ model generally applies to diseases where infection does not confer long-term immunity, such as the common cold, influenza, or herpes. Considering the repeated emergence of new viral variants of SARS-CoV-2 (BA.1, Omicron, EG.5...) that render existing immunity less effective, and knowing that previously infected people can become infected again, the $SIV$ model is a reasonable minimal model for the SARS-CoV-2 epidemics.

To study the effects of vaccination in the $SIV$ model, we further divide the population into a third fraction, $V$, representing the vaccinated people (Figure 1A). Two parameters associated with $V$ are $\bar{\beta}$, the breakthrough infection rate (the likelihood of a vaccinated individual becoming re-infected), and $\phi$, the willingness of a susceptible individual to take the vaccine. Considering the non-trivial feedback of the epidemic into the vaccine willingness, we assume that the probability of vaccination uptake, $\phi$, depends on the risk perception of infection and vaccination.

We use the simplest form of a vaccination game ^[21] under the presumption of uniform information exchange. Susceptible individuals can choose one of two strategies, either vaccinate ($V$) or not. The latter means they take a chance to become infected ($I$). The payoff of the strategies chosen by a susceptible individual depends on the relative strategies selected by others in the population. Denoting $r_{v}$ and $r_{i}$ as the risk perceptions associated with vaccination and infection, respectively, the equation for how $\phi$ changes can be expressed as

Here $\Delta E = -r_{v}V+r_{i}I$ is the marginal expectation difference between the two strategies ^[21,22,23]. Because scaling the payoff matrix does not change the qualitative dynamics of the game ^[23,24], the marginal expectation difference can be expressed as

where $f\equiv \frac{r_{i}}{r_{i}+r_{v}}$ is dimensionless. For simplicity, we call $\Delta E$ the "fear factor" and that plays a similar role to the incentive to vaccinate, while the fractions $f$ and $1-f$ are associated with the relative fear of infection and vaccination, respectively. We will show that $\Delta E$ has a significant impact on the oscillatory behavior of infection and vaccination.

The system of equations of the $SIV$ model is:

Because we do not consider immigration, population growth, and death, the population is conserved, $S+I+V = 1$, and the model can be reduced to three state variables by rewriting $S = 1 - (I + V).$ Figure 1 illustrates the box-and-arrow diagram and a triangular prism representing the three-dimensional (3D) phase space for the solutions.

The system of three ordinary differential equations, involving $I, V,$ and $\phi$, is non-trivial to analyze, even though it is straightforward to show that trajectories starting on the boundary (a side, edge, or vertex) of the prism can either remain in those domains or move into the prism. We study this model numerically and analytically to determine how $I$, $V$, and $\phi$ change over time within the prism. We limit ourselves to a reasonable estimate of the infection rate for susceptible, $\beta = 0.15$ (in units of 1/day), the infection recovery rate $\gamma = 0.07$ (1/day) for COVID-19, such that $R_0 = 2.1$ for the earlier variants of SARS-CoV-2 ^[25,26].

To study the effects of the vaccine, we consider a breakthrough infection and focus on the dynamics of the epidemic due to vaccine hesitancy by varying the fear factor $f$. Numerous studies estimated the breakthrough infection rate for SARS-CoV-2 vaccines ^[27,28]. In the US, UK, Brazil, South Africa, and Israel, the breakthrough infection rates have been estimated to be roughly 40% to 77% of the normal infection rate, hence our choice of breakthrough infection rate of $\bar{\beta} = 0.05$. It should be noted that the breakthrough infection depends on a multitude of factors, including vaccine type, number of doses, time after vaccination, virus variant, and individual immune response (older individuals have weakened immune systems and are more likely to become re-infected) ^[29,30,31]. Thus, varying $\bar{\beta}$ coincides with exploring different scenarios with variable vaccination efficacy due to the aforementioned factors. The incorporation of breakthrough infection is applied to consider the timeline for waning infection and reduced transmissibility between vaccinated and susceptible individuals. Indeed, the likelihood of transmission is reduced for vaccinated individuals ^[32], yet it could depend on the variant or the number of vaccines received, and age group ^[30,31]. Our choice of $\bar{\beta}$ lies between the intrinsic recovery rate $\gamma$ and the normal infection rate $\beta$. Unless otherwise stated, we will use these baseline parameter values (Table 1) in our numerical simulations.

3.   Three branches of asymptotic behavior

Seeking steady states, we set the right-hand sides of Eqs (2.1)–(2.4) to zero. These define the nullclines of 2D surfaces within the 3D prism. The most relevant nullclines to explore vaccine hesitancy are associated with $\dot{\phi} = 0$, namely, $\phi = 0, \; \phi = 1$, or $\Delta E = 0$. These nullclines correspond to three branches of distinct behaviors. We label the bottom plane of the prism ($\phi = 0$) as the "total rejection" branch and the top plane ($\phi = 1$) as the "full acceptance" branch. More generally, when the epidemic is perceived to be severe, there is less vaccination hesitancy, and $\phi$ increases. Conversely, if the epidemic is mild and few are infected, the perceived fear of infection is low, and $\phi$ decreases. Because $\phi$ is driven up or down by the sign of $\Delta E$, oscillatory behavior can quickly emerge from this feedback mechanism between the fear of infection and the fear of vaccination. Therefore, we label this third part of this nullcline ($\Delta E = 0$, a vertical plane in our prism) as the "oscillation" branch.

Along with $\dot{\phi} = 0$, the other nullclines ($\dot{I} = 0$ and $\dot{V} = 0$),

intersect at fixed points, with values denoted by $I^{\ast }, V^{\ast },$ and $\phi^{\ast }$ (and $S^{\ast } = 1-I^{\ast }-V^{\ast }$) for convenience. We will examine the stability properties of these fixed points by their respective Jacobian matrix:

3.1.   Total rejection: the "bottom" plane

This plane is characterized by $\phi \equiv 0$, i.e., no susceptible individuals are willing to be vaccinated. Thus, if no one in the population is vaccinated to begin with, $V$ will remain zero. However, even if the system is initiated with some fraction of $V$'s, that fraction can only decrease (via breakthrough infections) and end up at $V = 0$. With $V$ removed from the system, the dynamics reverts back to that of the trivial $SIV$ model, i.e., ending at the endemic stable steady state (provided $\beta > \gamma$):

We illustrate this scenario in Figure 2, where multiple initial conditions converge to this stable endemic state. Exceptions include the unrealistic case of 100% efficacy ($\bar{\beta} = 0$) and when $I$ exponentially decays with time $t$ and the entire population becomes vaccinated and permanently immune (equivalent to $R$ in a $SIR$ model). In the latter case, $I$ is so mild that the epidemic would end even without vaccination.

Next, we study small perturbations from this plane, i.e., $\phi \ll 1$. It is clear from Eq (2.4) that such perturbations will either decay or grow, depending on the sign of $\Delta E$. For example, the likelihood of a susceptible individual will choose to be vaccinated increases when there are many infected individuals ($I$), as the fear of infection becomes greater than the fear of vaccination. Quantitatively, we carry out a linear stability analysis and note that $\mathbb{J}$ is triangular with eigenvalues $-\beta I_{0}^{\ast }, -\beta I_{0}^{\ast }, fI_{0}^{\ast }$. Thus, although the bottom plane is unstable, this fixed point is stable within the bottom plane. Trajectories approaching the rest of the plane (away from the fixed point) will either be attracted or repelled to the plane, depending on the sign of the fear factor $\Delta E$.

3.2.   Full acceptance: the "top" plane

The total acceptance ($\phi \equiv 1$) branch corresponds to the top plane of the prism. Trajectories that start from this plane never leave and, as may be expected, always end at the $V = 1$ vertex. To confirm this scenario, we first consider the nullclines to seek a steady state in the plane.

The second part of Eq (3.1) specifies a straight line that joins the point $S = \gamma /\beta$ on the $V = 0$ edge and the point $V = \gamma /{\bar{\beta}}$ on the $S = 0$ edge. Because we have chosen $\beta, \gamma$ so that the former lies in the physical region, but the latter would lay outside if $\bar{\beta} < \gamma$ (indication of a good vaccine), this nullcline intersects the $I = 0$ edge at $S_{1}^{\ast } = \left(\gamma -\bar{\beta} \right) /\left(\beta -\bar{\beta}\right) = 1-V_{1}^{\ast }$. Meanwhile, since $I = 0$ is also a nullcline, the physically meaningful $\dot{I} = 0$ nullcline consists of both line segments: $I = 0$ and $I = I_{0}^{\ast }-\left(1-\bar{\beta}/\beta \right) V$ (red lines in Figure 3A), joined at $\left(S_{1}^{\ast }, V_{1}^{\ast }\right)$.

The second nullcline is a hyperbola, as $\left(\bar{\beta}I^{\ast }+1\right) \left(\bar{\beta}V^{\ast }+1\right) = 1+\bar{\beta}$ is an equivalent form for Eq (3.2). Symmetric in $IV$ and "anchored" at the $I$ and $V$ vertices of the triangle, it bows out to a mid-point $I_{m} = V_{m} = \left. \left[ \sqrt{1+\bar{\beta}}-1 \right] \right/ \bar{\beta}$ (blue line in Figure 3). For parameters of interest ($\bar{\beta} < \gamma < \beta$), the two nullclines meet at the stable steady state ($V = 1$ vertex). Figure 4 provides an example illustrating this plane's typical evolution trajectories.

Mathematically, the $V$ nullcline (blue hyperbola) can intersect the $I$ nullcline (red line) at two other points (Figure 3B). This scenario can occur only for $\beta > \left(1+\bar{\beta}\right) /2$. Keeping $\lambda = 0.07$ and $\bar{\beta} = 0.05$, the minimum value of $\gamma / \beta = 0.022$ for two intersections. This value corresponds to a $R_0 = 45.5$, much larger than any known reproduction number of SARS-CoV-2. This number corresponds to an extremely contagious infection, much higher than for measles ($R_0 > 10$ with a wide range ^[33]) or even the most contagious Omicron sub-variants ($R_0 = 9.5$ ^[34]). Because our main interest is to study a population with vaccination choices and how it affects "moderately" transmissible pathogens, we will not analyze these possibilities further.

Likewise, the two nullclines can intersect outside the feasible region if $\bar{\beta} > \gamma$. This scenario corresponds to an ineffective vaccine that cannot end the epidemic even if the entire population accepts it. We do not analyze these scenarios and summarize the obvious: with full participation ($\phi = 1$), an efficacious vaccine ($\bar{\beta} > \gamma$) will terminate an epidemic as long as the transmissibility is not too high ($\beta < \left(1+\bar{\beta}\right) /2$).

Within the top plane, every individual "accepts" the vaccine ($\phi = 1$) and the stable steady state is trivially the case of $V = 1$. Such "total acceptance" is best interpreted as a compulsory vaccine (such as polio). Clearly, for this population, the fear factor $\Delta E$ plays no role as everyone takes the vaccine regardless of social factors such as risk perception and fear of unwarranted side effects. Not surprisingly, trajectories on this plane that start with either sign of $\Delta E$ converge robustly to $V = 1$ (Figure 4). By contrast, $\Delta E$ plays a significant role away from this plane, as trajectories are attracted either towards the plane (when $\Delta E > 0$) or away from it (when $\Delta E < 0$).

3.3.   Oscillation: the "vertical" plane

The third and most interesting branch is where trajectories oscillate across the $\Delta E = 0$ surface (the green plane in Figure 1B). A steady state represents a perception of balanced fears that satisfies

This equation leads us readily to the balanced fears steady state (BFSS):

Note that these values are all physical (positive fractions) for the "realistic'' scenarios ($\beta \geq \gamma \geq \bar{\beta}$). In particular, $S^{\ast }\in \left[ S_{0}^{\ast }, S_{1}^{\ast }\right]$ as $f$ varies from $0$ to $1$.

To understand how the BFSS changes as the fear of infection increases from $0$ to $1$, we start with the $V = 0$ edge on the bottom plane: If someone has no fear of being infected, why would they bother taking a vaccine? Indeed, as $f \rightarrow 0$, this fixed point approaches the one in the "total rejection" plane, which is equivalent to the $SIV$ model (Table 2). On the other hand, when $f = 1$, the infection is eradicated through mass vaccination. This corresponds to timely vaccination interventions when the population is educated on the health benefits of vaccines and vaccine hesitancy is suppressed. Algebraically, it is surprising that this fixed point returns to this plane as $f\rightarrow 1$ (since $\phi ^{\ast }\propto 1-f$).

To appreciate this scenario, consider a population with a strong fear of infection. They are more likely to accept the vaccine, driving $V$ up and $I$ down. Eventually, the vaccination prevalence is high enough such that $\Delta E$ changes sign and $\phi$ decreases and reaches a small value of $\phi ^{\ast }$. If the population is in a regime where the proportion of vaccination exceeds the susceptible such that $\beta S+ \bar{\beta}V < \gamma$, then infections cannot regain their foothold, and the epidemic ends with $I\rightarrow 0$. Indeed, a more interesting behavior is associated with this BFSS as we study its stability properties next.

Consider the characteristic equation $\mu ^{3}-T\mu ^{2}+C\mu -D = 0$ for the matrix $\mathbb{J}$ with trace ($T$), cofactor ($C$), determinant ($D$), and the discriminant of the cubic equation $\Delta$:

We see both $T$ and $D$ are negative everywhere, except for $f \sim 0$. Thus, we encounter a pair of complex eigenvalues, which indicate oscillatory behavior. The physical interpretation is clear: with any imbalance ($\Delta E \neq 0$), $\phi$ is changed in order to restore that balance. But then, the system overshoots, leading to oscillations. The transition from a stable focus to an unstable one (a Hopf bifurcation) shows trajectories turn from spiraling into the fixed point to outward spirals ending in a stable limit cycle.

To analyze this bifurcation, let us label the eigenvalues of $\mathbb{J}$ by $\Lambda, \lambda \pm i\omega$. In terms of these, we have $T = \Lambda +2\lambda$, $C = 2\Lambda \lambda +\lambda ^{2}+\omega ^{2}$, and $D = \Lambda \left(\lambda ^{2}+\omega ^{2}\right)$, so that we can conclude $\Lambda < 0$ from $D < 0$. Further, for the regime where $\omega$ is real, the sign of $\lambda$ is conveniently identical to that of the characteristic $\Xi \equiv TC-D$. Thus, by solving $\Xi \left(\gamma, \beta, \bar{\beta}, f_{c}\right) = 0$, we can determine $f_{c}\left(\gamma, \beta, \bar{\beta}\right)$, i.e., the critical $f$ as a function of the parameter set $\left(\gamma, \beta, \bar{\beta}\right)$. For example, with our choice of $\left(\gamma, \beta, \bar{\beta}\right) = \left(0.07, 0.15, 0.05\right)$, we find $f_{c}\simeq 0.58$. For $f > f_{c}$, the trajectories in our model end either in limit cycles or, for more extreme values of $f$, in fixed points in the boundary of the prism. To illustrate this transition, we present two cases: $f = 0.50, 0.59$ (Figures 5 and 6), where the trajectory spirals into a stable fixed point and a stable limit cycle.

At the Hopf bifurcation, the frequency of the emerging limit cycle is $\omega = \sqrt{C}$, evaluated at this point. Though an explicit expression for $\omega \left(\gamma, \beta, \bar{\beta}\right)$ is available, we do not reproduce it here, as it is cumbersome and provides little insight. Its numerical evaluation agrees with the results of simulation runs. In particular, we find that the period of the limit cycle starts at about 120 days for $f\gtrsim f_{c}$, rising by a factor of about 3 when $f\simeq 0.9$.

Far above $f_{c}$, it is common to find cycles in which the system spends a long time near the $\phi = 0$ plane. At some point in the cycle, with substantial $V$ and little $I$, $V$ trickles out very slowly (as we have a low rate of breakthrough infection, $\bar{\beta}$) in a quasi-linear fashion so that $I$ increases slowly. As soon as the $I/V$ ratio passes $\left(1-f\right) /f$, the point of balanced fear, $\phi$ rises precipitously. This rise drives $V$ to increase rapidly (at the expense of $S$). As $V$ surpasses $I$, the population changes its perception of the relative risks when $(1-f)V > fI$ ($\Delta E$ flips sign) and $\phi$ swings downwards. This swing is also relatively fast as $\phi$ approaches zero. An example is provided in the inset of Figure 5A. From a game theoretical perspective, the changing of signs of $\Delta E$ corresponds to the emergence of the non-trivial Nash equilibrium, where a small perturbation about $\Delta E = 0$ causes an oscillation in the frequency of strategy selection. The essential properties of such swings can be understood through a solvable toy model, which demonstrates the nature of the $\phi$ equation (details in Appendix).

In a slightly more general form of the model, we can write $\dot{\phi} = \delta \phi\left(1-\phi\right) \Delta E$, where $\delta$ corresponds to the "strength of belief" and models how rapidly susceptible individuals get vaccinated when they decide to do so. The strength of belief mainly impacts the characteristic time of the oscillations (Figure 7B). The higher $\delta$ results in more rapid changes in $\dot\phi$ with oscillations in $V$ and $I$ in a more square shape.

In other highly infectious diseases with different $R_0$ values, as the fixed points in all three branches of the solutions depend on $\gamma / \beta = 1/R_0$, changing $R_0$ will shift the steady state values without changing the overall dynamics. Most notably is the critical fear factor $f_c$ for oscillation depends on $R_0$.

Even when between 60% and 90% of the population is vaccinated, because of the low vaccination efficacy ($\bar{\beta} \approx \gamma$), the breakthrough infection leads to fast spikes of infection followed by high vaccination uptake in response to the epidemic. This resembles the typical epidemic characteristic after a new viral variant emerges. Note that $\bar{\beta} > \gamma$ occurs only in the case of an ineffective vaccine.

4.   Summary and outlook

The COVID-19 pandemic caused over 7 million deaths worldwide and 1.2 million in the US alone (as of June 2024), more than vehicle crashes, gun violence, the flu, or many other public health threats combined ^[35,36]. The situation is unfortunate since many of these deaths could have been prevented if only more people had received vaccine shots. A recent study estimated that at least 232,000 deaths could have been prevented among unvaccinated adults during the 15 months in the US (from May 30, 2021 to September 3, 2022), had they been vaccinated with at least a primary series of vaccines ^[37]. Hence, accurate communication of valid scientific findings to the public and a decrease in the politicization of science are necessary.

Our minimal $SIV$ model adds a vaccine willingness to an $SIV$ epidemic model and explores how vaccine hesitancy impacts the epidemic evolution. A substantial body of literature exists for more complex models that account for specific COVID-19 characteristics and data. Unfortunately, these complex models are often difficult to analyze and understand in terms of the underlying complex dynamics. For example, at the onset of the COVID-19 pandemic, many epidemiological models were rapidly adapted to emerging public health data and provided valuable support for policymakers to develop public health responses. In particular, the Center for Disease Control and Prevention (CDC) continues to use a large and growing list of such models to provide forecasts of new cases, hospitalizations, and deaths (https://www.cdc.gov/coronavirus/2019-ncov/science/forecasting/mathematical-modeling.html).

Our simple model allowed for detailed analysis to shed new insights into how vaccine hesitancy can be a novel source of periodic oscillations in epidemic spread. We demonstrated that when vaccine willingness depends on the epidemic size, the disease-free equilibrium is only achievable in the case of total vaccine acceptance. On the other hand, in the case of total vaccine rejection, the epidemic invariably arrives at the endemic equilibrium with a finite infected population depending on the infection rate and recovery rate. In the more realistic case of intermediate vaccine hesitancy, the relative strength between the fear of infection and the fear of vaccine drives the epidemic progress. When the fear of infection is low, the epidemic settles into a steady state, with the number of people vaccinated proportional to the fear of infection and the number of people infected proportional to the fear of vaccination. Increasing the fear of infection, the system will cross a point where the fear of infection and vaccination balance out, and the system undergoes a Hopf bifurcation and enters periodic oscillations, eventually settling into a limit cycle. Many other studies have shown the existence of oscillations in COVID-19 modeling ^[38,39,40], but none of them studied vaccine hesitancy. Our minimal model offers novel insight into the role of vaccine hesitancy as a novel source of epidemic oscillation, and the challenges of achieving herd immunity via vaccination in the presence of vaccine hesitancy.

There are innumerable avenues worthy of future research; only a few will be mentioned here. The present model has minimal parameters associated with various transition rates. Other rates that may affect the qualitative outcomes of the model include, e.g., births/deaths and waning immunity. In the latter, $V\rightarrow S$ with an additional timescale. Similarly, we may incorporate more sophisticated ingredients to model the behavior of vaccine acceptance. Another extension of our simple model would be to include stochastic effects or formulate it as an agent-based model to incorporate an individual's history. Preliminary studies involving agent-based simulations revealed several interesting results, e.g., finite lifetimes of the endemic steady states (also known as lifetimes to extinction in the language of population dynamics).

Our model consists of individuals who have the same information and use the information to assess relative risks, which is a simplification of information processing, selective exposure, and intrinsic bias due to one's political affiliation, media exposure, and social network ^[41,42,43]. A separate, network-based model for information transmission in addition to a physical network-based disease transmission would be desirable to capture the more realistic information dynamics.

It is also possible to study fluctuations of and correlations between modeling parameters (e.g., transmission rates, recovery rate, and fear factor); both are absent from the analysis in this paper. Most diseases, especially SARS-CoV-2, affect different age groups differently. We can investigate in more details how an epidemic might spread and end by including age distributions. In general, only some people are in contact with everyone else in the population due to social and geographic constraints. Overlaying social and spatial structures on our model and investigating their effects could provide new insights. Finally, using agent-based models, we consider the much more realistic situation where there is a variety of responses (individual $\phi$'s) due to different levels of fear and misinformation.

We believe that understanding how vaccine hesitancy and reluctance impact the progression of an epidemic can help to design better vaccine campaign strategies and shorten the duration of an epidemic.

5.   Appendix

An analytically solvable toy model

Here, we show a minimal model that can be solved analytically and may provide some insight into the oscillatory behavior in our model. Consider Eq (2.4) and the simplified equation for $\Delta E$:

where $A, B, \sigma$ are positive constants. Although these constants can be absorbed into scale factors of the variables and $t$, their explicit presence facilitates discussions here.

This system of equations can be solved exactly, as in classical mechanics. To make the analogy explicit, let us think of a particle with unit mass moving in $x\in \left[ -\infty, \infty \right]$. In addition, motivated by $dx = d\phi /\phi \left(1-\phi \right)$, let us change variables to

(An arbitrary $x_{0}$ may be added to $x$, we choose it to be $0$ when $\phi = 1/2$.) Then

lies in $\left[ 0, 1\right]$. Equations (5.1) and (5.2) become

which displays classic Hamiltonian dynamics with

and the potential

For $\left\vert x\right\vert \gg x_{0}$, the potential $\mathcal{V}$ is almost linear, being $-\sigma Bx$ for $x < 0$ and $\sigma Ax$ for $x > 0$. Keeping $A, B$ constant and varying $\sigma$, we see that $\sigma$ corresponds to the steepness of the V-shaped potential. We can verify that, given the same initial conditions ($\Delta E\left(0\right), \phi \left(0\right)$), a large $\sigma$ will result in $\phi$ making rapid transitions from $0$ to $1$ (i.e., our particle oscillating between large positive and negative values of $x$). Meanwhile, the "momentum" changes sign at the bottom of $\mathcal{V}$ almost instantaneously, corresponding to the kink-like structures at extreme values of $\Delta E$. When the trajectories take us to large values of $x$, the almost linear potential $\mathcal{V}$ resembles uniform gravity, leading to the momentum being linear in $t$. This aspect captures the essence of quasi-linear regions in a $V-I$ plot.

Finally, the main advantage of this toy model is that it guarantees limit cycles. Indeed, this is a neutral system with an infinite number of limit cycles, labeled by the value of $\mathcal{H}$, the total "energy" of our particle.

Steady states

We compute the steady states by setting Eqs (2.1)–(2.4) to zero and solving for $(S^*, I^*, V^*, \phi^*)$. Our model demonstrates that the $SIV$ model can be reduced to the $SIV$ framework in the extreme cases of zero or all vaccinations. We also include a unique complex balanced fear steady state (BFSS). The steady states can be illustrated in the phase portrait (Figure 8) and the corresponding Table 2. Note the steady states and their stability depend on the initial conditions and parameter values.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was partially supported by Public Health Service grants 1R01EY028450 from NIH/NEI, the endowment for the Frady Whipple Chair in Mathematics at Georgia State University (to Y.J.), the endowment for the Evelyn and John G. Phillips Distinguished Chair in Mathematics at Tulane University and R01HD086794 from the NIH/NICHD and U01GM09658 from the NIH/NIGMS/MIDAS (to JMH), and a 2CI PhD Fellowship from Georgia State University (to AM). RKPZ acknowledges valuable discussions with Kevin Bassler and Hedda Grelz.

Conflict of interest

The authors declare there is no conflict of interest. Yi Jiang is an editorial board member for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article.