Incorporating changeable attitudes toward vaccination into compartment models for infectious diseases

1.   Introduction

Formulating public health policies toward infectious diseases requires the capability to predict how the infectious population will change over time under various scenarios. Accurate modeling of the effects of public health measures is essential. During the epidemic phase of a novel infectious disease, public health measures typically include basic precautions such as quarantine, masking, and isolation. These tools become less important after a disease enters the endemic phase, where a large part of the population is no longer susceptible and infectious populations are a small fraction of the total population. From then on, vaccination becomes the most important public health tool, particularly for diseases where immunity wanes over time, provided a vaccine is available.

Compliance with public health measures is not compulsory in most societies; hence, noncompliance can be a significant factor that reduces the effectiveness of those measures. In 2013, Dubé et al. ^[1] used sociological data to study this phenomenon in regard to the increased hesitancy at that time among parents toward vaccination for the standard childhood diseases. The authors quantified individual attitudes toward vaccination on a continuum from outright refusal to full acceptance. An individual's location within this continuum is informed primarily by their perception of the risk of the disease (and the vaccine) and the perceived importance of vaccination in mitigating that risk. These perceptions are based on the individual's background knowledge, acquired information (accurate or not), and cultural norms. Acquired information can come from public health announcements, recommendations of health professionals, communication with acquaintances, journalists, and social media posts.

Dubé's findings indicate a need for mathematical modelers to incorporate opinion dynamics into epidemiology models. An early paper on this topic examines the impact of human behavior on a disease outbreak for which control is limited to the basic mitigation procedures of quarantine and isolation ^[2]. This paper does not consider these behaviors to be subject to hesitancy, but assumes that the primary driver of opinion is awareness of the disease outbreak. The model categorizes individuals into groups based on different awareness levels. An individual's awareness level can increase through information transmission via personal communication and fade over time through a natural decay process.

Other approaches to combined epidemiological and opinion dynamics have appeared since the COVID-19 pandemic highlighted the extent of vaccine noncompliance. In recounting these developments, we consider only those papers whose models include vaccination. Ali et al. ^[3] use an epidemiological model for smallpox as part of a bioterrorism scenario. The vaccination component categorizes people as cooperative or non-cooperative. Opinion dynamics are driven by encounters between individuals, with no consideration of public health or social media communication, and are assumed to result in increased cooperation. A related paper ^[4] couples a COVID-19 model with an opinion dynamics model similar to that of Ali et al. ^[3], but with an additional mechanism for the influence of public communications. The authors assume that interpersonal and social communications increase cooperation.

Another trend is to use an opinion dynamics model based on statistical physics; this formulation treats opinion as a continuous variable that changes through random interactions with other individuals and stochastic variation ^[5,6,7,8,9]. With no processes involving public communication or epidemiological conditions, the result is a model that is mechanistically correct for the random drift process but neglects processes that are likely to be much more important than random drift. Model outcomes are largely determined by the initial distribution of opinions.

The papers cited above do not include a mechanism for opinions to respond to the epidemiological state, which we consider the most quantifiable driver of opinion change. A seminal paper by Bauch ^[10] does include such a mechanism by assuming that the fraction of individuals who get vaccinated in a SIR model where vaccination confers immunity is subject to a dynamic game theoretic model based on the idea that individuals determine whether to accept vaccination by making a rational choice comparing the perceived benefits of vaccination and non-vaccination. Individuals change their status when information, tracked using an information index, is exchanged with other individuals, with an overall increase in vaccination behavior when disease prevalence exceeds a certain threshold. d'Onofrio et al. ^[11] extend Bauch's model with more general functions for the status change dynamics, as well as functions with a time delay. In a similar model ^[12], the same authors add public information as a contributor to the information index.

Models developed by d'Onofrio et al. ^[13] and Buonomo et al. ^[14] also use an information index, but have it directly impacting the vaccination rate rather than changing rates through individuals interacting. d'Onofrio et al. incorporate the information index into a SIR model and Buonomo et al. do so in an SEIR model, both assuming that vaccination confers lifelong immunity.

The onset of COVID-19 has influenced the direction of mathematical epidemiology in numerous ways, including the heightened importance of individual attitudes toward vaccination and the phenomena of imperfect vaccination, in which a vaccine decreases the rate of transmission without conferring immunity. Buonomo et al. ^[15] incorporate an information index model into a COVID-19 epidemiological model with imperfect vaccination. The information index changes through a process of linear decay toward an equilibrium proportional to the infectious population, and the vaccination rate increases monotonically with the current value of the index. The model is shown to have a unique endemic disease equilibrium when the disease-free equilibrium is unstable, but general stability results are not established. In another recent paper, Zuo et al. ^[16] return to the game-theoretic approach of Bauch et al. ^[10] and d'Onofrio et al. ^[11], but with a base epidemiological model that assumes imperfect vaccination. There is no existence or stability analysis for endemic disease equilibria performed in that paper.

Xuan et al. ^[17] also feature vaccination behavior dependent on disease prevalence. The authors use an agent-based SIS model with each individual's avoidance behavior based on their opinion attribute. Changes in this attribute over time are based in part on infection probability and in part on a process of consensus-building among agents connected in the social network. Because the disease state is such an obvious driver of attitude toward vaccination, it is worthwhile to consider models in which this single driver of opinion dynamics is added to a relatively simple epidemiological model.

All of the papers cited above incorporate vaccination into the model as a single-phase spontaneous transition from a susceptible class to a recovered class ^[18,19,20]. This type of model implicitly assumes that everyone will be vaccinated unless they are infected first and that there are no limits on vaccine supply and distribution. In practice, however, a significant fraction of people do not get vaccinated for a variety of reasons. Some are unable to be vaccinated because of other health conditions, some are committed opponents of vaccination, some are hesitant because of concerns about safety, efficacy, or necessity, and others are just insufficiently motivated. This phenomena is sufficiently common that any model that incorporates vaccination ought to account for it.

Ledder ^[21] addresses these issues by proposing realistic models for vaccine supply and distribution in the epidemic phase and by partitioning the susceptible class into Pre-vaccinated and Unprotected subclasses, the former consisting of people who are willing to be vaccinated but are waiting for vaccine availability, and the latter for people who are either unwilling or unable to be vaccinated. In the epidemic phase, when vaccines may not yet be readily available, the presence of a subclass of the unwilling and unable has only minimal impact due to the existence of a pool of willing individuals awaiting vaccination. However, the existence of an unprotected subclass is an important driver of the long-term disease dynamics.

Ledder's models ^[21] assume that the fraction of unwilling/unable individuals is fixed over time. In reality, attitudes towards vaccination and other public health measures may change in response to varying circumstances. In this paper, we aim to develop a model that addresses vaccine hesitancy by partitioning the population into Willing and Unwilling/Unable subclasses, similar to the approach in Ledder ^[21]; however, we also include processes whereby individuals can move from one of these subclasses to the other due to a change in opinion.

Our goal in this paper is not to create a realistic model for any particular disease, but to create a model that incorporates the most essential features of disease dynamics along with vaccination, gradual loss of immunity, and subdivision into Willing and Unwilling/Unable subclasses with transition dynamics based on disease prevalence. To that end, we begin with an SIR model in which each of the three classes is divided into Willing and Unwilling/Unable subclasses. We allow the parameters that control the attitude-change processes to be nonlinear functions of disease prevalence, giving specific forms for these functions only when needed for examples. We will see that the nonlinearity in the functions that determine the attitude-change rate parameters can be sufficient to yield instability.

For simplicity, we assume perfect vaccination in the sense that vaccination confers the same degree of immunity as recovery and that individuals in the Willing Recovered subclass receive vaccination updates sufficient to maintain immunity. In addition to omitting more subtle vaccination effects such as partial effectiveness or partial immunity, we omit disease-induced mortality. While these features would broaden the applicability of the results, they would also make the model more complicated without changing the qualitative behavior.

We begin in Section 2 with the development of the model and its reformulation with dimensionless variables. Section 3 presents the analytical results for the model. Section 4 contains numerical simulations to illustrate the results and a global sensitivity analysis showing the effect of relevant parameters on long-term prevalence. Finally, we conclude with a discussion of the results and their practical implications in Section 5.

2.   Model development

We begin by presenting an initial model, which we then scale to provide a more convenient version for analysis. The model also requires two functions to describe how the stance towards vaccination depends on the infectious population; some example functions are given at the end of this section. Notation for the model is summarized in Tables 1 and 2. Note that all parameters must be nonnegative for the model to represent an epidemiological scenario.

2.1.   The initial model

Our model extends the PUIRU model presented by Ledder ^[21] by incorporating additional processes that allow individuals to change their status from Willing to Unwilling/Unable and vice versa. Figure 1 shows a schematic diagram of our model, with definitions of relevant variables and parameters in Table 1. The model is based on the following assumptions:

1) The population can be divided into compartments representing the Susceptible ($S$), Infectious ($I$), and Recovered ($R$) population fractions. Each is divided into two subgroups: a 'Willing' subclass of people who get a vaccination when it is recommended (subscript $W$), and a complementary 'Unwilling/Unable' subclass for those who are either unwilling or unable to be vaccinated (subscript $U$).

2) Susceptible individuals become infectious through contact with infectious individuals with rate constant $\beta$.

3) Infectious individuals recover through a single-phase transition process with rate constant $\gamma$.

4) Individuals in the Willing Susceptible subclass ($S_W$) move directly to the corresponding Recovered subclass ($R_W$) upon vaccination, with rate constant $\phi$.

5) Individuals can change from Willing to Unwilling/Unable subclasses and vice versa with rate coefficients that are functions of the total Infectious population: $\Psi(I)$ and $\Omega(I)$, respectively.

6) Unwilling/Unable Recovered individuals ($R_U$) lose immunity over time, with rate constant $\theta$.

7) Willing Recovered individuals ($R_W$) remain in that subclass because they receive booster vaccination doses as needed, except for those who switch from Willing to Unwilling/Unable at rate $\Psi(I)R_W$.

8) All newborns are assumed to be Unwilling/Unable Susceptible.

9) There is no mortality caused by the disease, either by direct fatality or indirect reduction of life expectancy. The natural birth and death rates are assumed to be equal, and the initial population size is set to ensure that the total population remains constant at a normalized value of 1.

We call particular attention to Assumption 5. By making the rate coefficients functions of the system's state, we introduce additional nonlinearity, which could produce features not observed in standard epidemic models where nonlinearity is limited to the infection process. This assumption adds a degree of additional realism to the model. The patterns of behavior seen in response to COVID-19 show that a lessening of the disease burden induces people to allow their immunity to lapse.

Note that individuals who change from Unwilling Recovered to Willing Recovered could in theory lose immunity before they receive their first vaccination. We neglect this possibility because vaccination occurs on a time scale of days to weeks, while loss of immunity occurs on a much slower time scale of months to years.

While it seems natural to write differential equations for each of the compartments in the diagram, it is more convenient to replace the differential equations for $S_U$, $I_U$, $R_W$, and $R_U$ with equations for $S$, $I$, $R$, and $W$ (the total population fraction in the Willing subclasses), resulting in the model

where $S_U = S-S_W$, $I_U = I-I_W$, and $R_U = (1-W)-(S_U+I_U)$. Note that the population size is conserved with a total population of 1; hence, the population fraction of Unwilling/Unable individuals appears in the $W$ equation as $1-W$.

2.2.   The reformulated model

Stability analysis is easier when models are scaled appropriately. The dependent variables are already scaled as fractions of the total population. We introduce a dimensionless time $T$ using the population time scale $1/\mu$, that is $T = \mu t$; hence,

We also introduce dimensionless parameters

and use $'$ to denote the derivative with respect to $T$. The intuitive meanings of all newly-defined notations here are included in Table 2. The parameter definitions deserve a comment. The parameter $\epsilon$ is the ratio of the demographic rate constant $\mu$ to the disease rate constant $\gamma+\mu$; hence, this parameter is very small and can be used for asymptotic approximation. The infection and vaccination processes are assumed to occur on the disease time scale; hence the rate constants $\beta$ and $\phi$ are scaled by $\gamma+\mu$. The loss of immunity and change of attitude processes are assumed to occur on the demographic scale, so their rate constants are scaled by $\mu$.

Note that the $R$ equation can be eliminated using $S+I+R = 1$, and therefore it is decoupled from the system. We reorder the equations to facilitate the use of the next generation matrix method for the calculation of the basic reproduction number in the presence of vaccination, $\mathcal{R}_v$. With these changes, the model becomes

The system (2.3) incorporates the naive assumption that all dependent variables are $O(1)$ quantities on the time scale of long-term population dynamics. The resulting system is not correctly scaled for analysis of long-term behavior, however, as the $S_W$ equation would then imply $S_W = O(\epsilon)$ at equilibrium, contradicting the inherent assumption that all variables are $O(1)$. With $S_W = O(\epsilon)$, the $I_W$ equation and $S$ equation then show that $I_W$ and $I$ are also $O(\epsilon)$. The solution is to build in the assumption that these three populations are in fact $O(\epsilon)$ on this time scale ^[22]. As noted above, $\epsilon$ represents the ratio of expected time spent in the infectious compartment to mean lifetime. Assuming infections last on the order of weeks and a lifetime is roughly 80 years, this parameter is indeed small. We therefore replace $I$, $I_W$, and $S_W$ by

At the same time, we introduce some lumped quantities that are convenient now or will become so later:

1) As needed, a bar over the top of a quantity indicates that quantity plus 1; for example, $\bar{\Sigma} = \Sigma+1$.

2) To simplify the equilibrium expressions, we remove a factor of $\mathcal{R}_0$ from $Y$, $Z$, and $P$ and redefine these variables accordingly:

3) Various sums of parameters arise in the analysis:

With these changes, we obtain the final form of the system:

Because the equations are now properly scaled, we can identify $O(\epsilon)$ terms as being relatively small compared to $O(1)$ terms. These $O(\epsilon)$ terms can be considered as small perturbations for a regular perturbation analysis, but here we consider only the leading order terms. This decouples the $Z$ equation in the system (2.7), leaving a four-component system that will be used for analysis:

2.3.   The status-change functions $\omega$ and $\psi$

Regardless of whether the scaled Unwilling/Unable-to-Willing transition rate coefficient $\omega$ is increasing or decreasing, it seems reasonable that the scaled Willing-to-Unwilling/Unable transition rate coefficient $\psi$ should be moving in the opposite direction. While one of these functions might be more sensitive to the infectious population fraction than the other, there is no way to know which. As a preliminary approximation, it is reasonable to assume that the sum $\Sigma = \omega+\psi$ remains approximately constant. To keep the model as simple as possible for the elucidation of general properties, we will think of $\Sigma$ as constant whenever we are ready to consider specific examples, specifying $\omega$ and replacing $\psi$ by $\Sigma-\omega$.

We can identify two different classes of plausible status-change functions $\omega$, given the assumption of constant $\Sigma$. When infection rates are low, $\omega$ should be increasing with $Y$. This is something that typically occurs with diseases; a lower incidence leads to less strenuous vaccination campaigns, as is currently done with some vaccinations administered only to people traveling to a country with known cases. As the infection rate increases, two types of behavior seem possible. Healthcare authorities will undoubtedly increase the emphasis they put on the need for vaccination, and an informed public that is trusting of healthcare authorities will be more likely to be vaccinated. On the other hand, public skepticism could lead to a situation where some people interpret higher disease burdens as vaccine failures and become less likely to be vaccinated. Thus, we consider functions that are monotone-increasing and functions that increase to a local maximum before declining.

To study general behavior, the specific families of functions used for the monotone and non-monotone cases can be chosen for mathematical convenience. In the sequel, we will use

as a family of monotonic functions and

as a family of non-monotonic functions. The latter family has the simple property that it vanishes to 0 as $Y \to \infty$; this seems reasonable, as having nearly everyone infected would likely make everyone lose confidence in the vaccine.

2.4.   Choice of parameter values

Values for parameters are problematic when a model is intended to be general rather than to match a specific disease. It is best to defer choosing values until absolutely necessary. In this case, we can obtain analytical results using only the approximation $\epsilon \to 0$, but parameter values are still required for visualization of results.

All of the dimensionless parameters depend on the principal fast and slow time scales. The value $1/\mu$ is the mean lifespan of the population, which we can safely take to be on the order of 80 years. The value $1/(\gamma+\mu) \approx 1/\gamma$ is the mean disease duration. While this can vary among diseases, it is reasonable to take 2 weeks as a typical value. With these choices, we have $\epsilon = 1/2000$, which means that results for realistic $\epsilon$ should be well approximated by results for $\epsilon \to 0$.

We consider a reasonable range of values for $\mathcal{R}_0$ to be from 2 to 10, with 10 a rough estimate for measles, and possibly also for current variants of COVID-19.¹ We use the intermediate value $\mathcal{R}_0 = 4$ for simulations.

¹The original strain had $\mathcal{R}_0 \approx 5.7$ ^[23], and more recent variants are surely higher.

The parameter $h$ represents the expected number of times immunity is lost during a normal lifespan. Actual values of this parameter can range from 0 for diseases like measles, where immunity is generally permanent, to higher values for diseases like the flu, where immunity wanes quickly. We will generally use $h = 0$ and $h = 10$ (mean immunity of 8 years) as representative of diseases with long and short immunity durations, respectively. In practice, larger values of $h$ could be justified for diseases like COVID-19.

While the parameter $\Sigma = \omega+\psi$ is variable in the general case, we consider it reasonable to treat it in practice as a fixed measure of community opinion flexibility. Individuals in a population vary in their willingness to reconsider their opinions, but we can take $\Sigma$ as a mean rate in the same way that other parameters are mean rates for incubation or disease duration. We use $\Sigma = 20$ as a practical upper bound; this value represents a population in which opinions change approximately 20 times per lifespan, or about every 4 years.

The parameter $v$ is the ratio of time spent in the infectious compartment to the time required for newly willing people to be vaccinated. The value $v = 0.5$ corresponds to a 4-week delay. The parameter $a$ from the nonlinear functions for $\omega$ indicates the sensitivity of opinion to disease prevalence. There is no data for this empirical parameter, but we will see that a full range of interesting results is possible without taking a value larger than 5. The parameter $d$ from the non-monotonic function for $\omega$ represents how close $\omega$ comes to its maximum $\Sigma$. A value of 0.9 serves as a reasonable upper bound.

3.   Analysis and results

We begin with the disease-free equilibrium (DFE), continue with general results for the endemic disease equilibrium (EDE), and then examine results for some specific attitude-change functions $\psi$ and $\omega$. Because this model is intended to study general features of epidemics, rather than reproducing data or making predictions for a specific disease, our primary interest is in qualitative results rather than quantitative details. For this reason, it makes sense to make use of $\epsilon \to 0$ asymptotics for the endemic disease equilibrium analysis.² Hence, terms in a sum that are $O(\epsilon)$ compared to the leading order term will be discarded. New $O(\epsilon)$ terms will also emerge in the stability analysis, where we will discard those that are directly comparable to higher-order terms.

²Analysis of the disease-free equilibrium does not require asymptotic approximation.

3.1.   The disease-free equilibrium and vaccine-reduced reproduction number

Computation of the vaccine-reduced reproduction number and identification of the stability requirement for the disease-free equilibrium are accomplished by standard methods. We present the result here and give details in Appendix A.

Proposition 1. The vaccine-reduced reproduction number for the model (2.8) is

The disease-free equilibrium is stable if and only if $\mathcal{R}_v < 1$; equivalently, there is a minimum value of $\omega$ needed for stability:

where $\rho(Y) = \bar{\Sigma} (Y)+ \mathcal{R}_0 h$.

Figure 2 illustrates the stability requirement for the DFE. Larger $\psi$ requires larger $\omega$ to achieve stability, and larger $h$ increases the $\omega$ requirement. As a practical matter, we should expect that this stability requirement cannot be achieved, as this would require a significant effort to vaccinate young children in spite of having no active disease in the population.

3.2.   The endemic-disease equilibrium

The parameters $\omega$ and $\psi$, and the quantities derived from them, are functions of $Y$. Thus, when we solve the equilibrium system for $y$, the implicit result that we get will be an equation that contains unknown values for $\omega$ and $\psi$ as well as $y$. This equation can be solved for $y$ only after the functions $\omega(Y)$ and $\psi(Y)$ have been defined. Given this context and the frequent occurrence of the combination $\Sigma = \omega+\psi$, it is practical to consider the equilibrium relation as defining a "target" $\omega$ value, which must be met at a specific positive value of $y$ to maintain equilibrium, considering the known values of $\Sigma$, $\mathcal{R}_0$, and $h$.

In the asymptotic limit $\epsilon \to 0$, with $Y > 0$, the system (2.8) yields an equilibrium system in which the $Y$, $W$, and $P$ equations reduce to

Substituting for $W$ and $p$ in terms of $\omega$, the $S$ equation becomes

which we can recast in terms of $\omega$ as

This result serves to provide a condition for finding equilibria---the value of $\omega$ prescribed by the model must match the value of $\omega$ that yields an equilibrium.

Proposition 2. For any given values of the parameters $\Sigma$, $\mathcal{R}_0$, $v$, and $h$, endemic disease equilibria occur when

where

with $r = \mathcal{R}_0-1$, $\bar{h} = h+1$, and $\rho = \bar{\Sigma}+ \mathcal{R}_0 h$.

Proof. Note that endemic disease equilibria must have $y < r \bar{h}$, or

The graph of $\omega^*$ is decreasing and reaches 0 at finite $Y$. Thus, monotone increasing functions for $\omega$ will have a unique EDE if and only if the disease-free equilibrium is unstable, that is, $\omega(0) < \omega_{cr}$. □

3.2.1.   Endemic disease equilibrium stability

While it is not normally possible to obtain analytical stability criteria for a 5-component system, such as (2.7), with unspecified parameter values, the 4-component system obtained with asymptotic approximation can be fully analyzed in the asymptotic limit (see Ledder ^[22] for a description of the method), which simplifies the criteria in addition to reducing the problem to four components. We summarize the results here, leaving details for Appendix B.

Proposition 3. Let $Y > 0$ be an equilibrium point for the system. Given the values of $\omega$, $\psi$, and their derivatives at the equilibrium point, the stability requirements for the endemic disease equilibrium in the asymptotic limit $\epsilon \to 0$ reduce to

where

For the simplified case where $\Sigma$ is independent of $Y$, these conditions reduce to

where

Figure 3 illustrates the stability requirements for $\omega'$ in terms of $Y$, with $\mathcal{R}_0 = 4$, $v = 0.5$, $h = 0, 2, 10$ (left to right) and $\Sigma = 20$ (solid) and $\Sigma = 5$ (dashed). The constant $\omega$ case corresponds to a model with no dependence of rates on disease level, for which the (unique by Proposition 2) EDE will always be stable. Instability results from greater sensitivity of $\omega$ to $Y$.

3.3.   Results for some specific status-change functions

We now consider results for some specific status change functions, all assuming $\Sigma$ is constant.

3.3.1.   Monotone functions

We begin with the broad class of monotone increasing functions for $\omega(Y)$. If $\omega(0) > \omega_{cr}$, then the equilibrium condition $\omega(Y) = \omega^*(\mathcal{R}_0 Y)$ will never be achieved, and the disease-free equilibrium will be stable. Instead, we assume $\omega(0) < \omega_{cr}$. In this case, the combination of properties $\omega' > 0$, $\omega(0) < \omega_{cr}$, and $\omega^*(r \bar{h}) = 0$ guarantees that there will be a unique endemic disease equilibrium $Y > 0$, and that it will be stable if and only if

with $\Upsilon$ given in (3.11).

As a simple example of this class, we consider the one-parameter family given in (2.9). Several functions in this family are shown in Figure 4, along with the curve $\omega^*(\mathcal{R}_0 Y)$ with $\mathcal{R}_0 = 4$, $v = 0.5$, $\Sigma = 10$, and either $h = 0$ or $h = 10$. The unique equilibrium for each curve is marked according to its stability result.

A broader summary of the $h = 0$ and $h = 10$ cases, with $\mathcal{R}_0 = 4$ and $v = 0.5$, appears in Figure 5, which shows level curves of $Y$ in the $\Sigma a$-plane along with the stability boundary. The stability boundary curves move to the left as $a$ increases; thus, larger $a$ is always destabilizing.

3.3.2.   Functions with a local maximum

We now turn to the broad class of functions for $\omega(Y)$ that increase to a maximum at a point ${Y_m < Y_{max}}$, and then decrease as $Y$ increases further.

When a specific example is required for illustration, we consider the family (2.10), which increases to a maximum value of $\Sigma d$ occurring at $aY = 1$. The requirement $\psi \ge 0$ translates into a parameter restriction $d \le 1$. A sampling of this family is illustrated in Figure 6. The dashed curve is the equilibrium condition $\omega(Y) = \omega^*(\mathcal{R}_0 Y)$, with $\Sigma = 10$, $h = 10$, and $\mathcal{R}_0 = 4$. All curves for $\omega$ have $a = 0.62$. The three curves have $d$ values of 0.5, 0.7, and 0.9, from bottom to top. The first of these has one (stable) equilibrium, with a value of $Y$ that is close to $Y_{max}$, where $\omega^* (\mathcal{R}_0 Y)$ is 0. The last has one (unstable) equilibrium, but this time with a small value of $Y$. The intermediate curve has three equilibria. As $d$ increases from 0.7 to 0.9, the second and third equilibria converge to a point and then disappear. Similarly, as $d$ is decreased from 0.7, the first and second equilibria converge and then disappear.

In general, suppose $\omega(Y)$ increases from an initial value $\omega(0) < \omega_{cr}$ to a maximum at $Y_m$ and subsequently decreases, maintaining $\omega(Y^*) > 0$ throughout. Then the equilibrium condition ${\omega(Y) = \omega^*(\mathcal{R}_0 Y)}$ will necessarily be achieved at least once, but may be achieved multiple times, with bifurcation points separating curves in the family that have different numbers of equilibria. The $\omega$ functions that mark bifurcations satisfy the requirement that the $\omega(Y)$ and $\omega^*(\mathcal{R}_0 Y)$ curves are tangent at a point, that is

Because ${\omega^*}' < 0$, bifurcation points can only occur in the domain where $\omega$ is decreasing, that is, when $Y > Y_{m}$; hence, each bifurcation point in the $ad$-plane corresponds to a unique $Y$ value in the interval $(Y_m, Y_{max})$, where $Y_{max}$ is defined in (3.5). The function $\omega_2$ in (2.10) has two bifurcation points for any given value of $a$. For example, for $a = 0.62$, the bifurcations occur at $d \approx 0.5507$ and $d \approx 0.7923$.

Comparison with earlier results discloses an interesting synergy, which holds for any function with the properties given above:

Proposition 4. Suppose $\omega(Y)$ increases to a maximum at a value $Y_m$, and then decreases, with ${\omega(0) < \omega_{cr}}$ and $\omega(Y_{max}) > 0$. Then for the region of solutions $Y > Y_m$, the stability boundary (3.10) coincides with the bifurcation curve (3.13).

Proof. The proposition follows immediately from the fact that $\omega' < 0 < \mathcal{R}_0 \Upsilon(\mathcal{R}_0 Y)$ when $Y > Y_m$, thus automatically satisfying the upper bound criterion for stability. □

Figure 7 illustrates the behavior of the system. Panel a shows the bifurcation curves (black), along with the stability boundary curve for positive $\omega'$ (blue), given by

For points above this stability curve, the smallest $Y$ value (the unique $Y$ outside of the region of multiple solutions) is unstable. When there are multiple solutions, the one with largest $Y$ is stable and the one with intermediate $Y$ is unstable, from Proposition 4. The regions are labeled with their stability pattern from low $Y$ to high, with the slender unlabeled region in the middle as SUS, that is, three equilibria that are stable, unstable, and stable, from smallest to largest. Panel b shows the solution curve in the $dY$-plane corresponding to $a = 0.64$. The $d$ axis is partitioned into four regions corresponding to the four stability patterns that occur for $a$ values in region 3 of panel a.

3.4.   Mean prevalence

In judging the public health impact of a disease, it is the long-term average prevalence $\bar{Y}$ that matters when the equilibrium solution is unstable. For parameter sets leading to an unstable EDE, $\bar{Y}$ can be determined numerically by augmenting the system of differential equations with a variable $Q$ that represents accumulated infectiousness, that is, it is governed by the differential equation

Now let $t_{j-1}$ and $t_j$ be the location of any two consecutive peaks in $Y$. The mean of $Y$ over the interval $[t_{j-1}, t_j]$ is

Assuming the solution is converging to a limit cycle, this quantity $\bar{Y}_j$ converges to $\bar{Y}$ as $j \to \infty$.

Figure 8 shows the results of some experiments for the monotone example $\omega_1$. Data for $\bar{Y}$ was computed by identifying the sequence of values $\bar{Y}_j$ and running the simulation until the differences among several consecutive values of the sequence are small enough to approximate convergence. Each panel compares $\bar{Y}$ to $Y^*$ over a range of values of one of the parameters $\Sigma$, $a$, $\mathcal{R}_0$, and $h$, using default values $\Sigma = 6$, $a = 2$, $\mathcal{R}_0 = 4$, $h = 10$, $v = 0.5$, and $\epsilon = 5\times 10^{-4}$. For oscillatory solutions, $\bar{Y}$ is larger than $Y^*$ for all examples we have studied.

The plots show an interesting difference in the behavior of the mean prevalence $\bar{Y}$ and equilibrium prevalence $Y^*$ with respect to the parameters. While $Y^*$ is monotone for each parameter (decreasing for $\Sigma$ and $a$, increasing for $\mathcal{R}_0$ and $h$), $\bar{Y}$ is not always so. Indeed, there is sometimes a local minimum with respect to one or more parameters that occurs at the stability boundary. It is particularly curious that it is possible for the mean prevalence to decrease with increasing infectivity $\mathcal{R}_0$.

A similar analysis cannot be done for the non-monotone case. Here, there are sometimes multiple stable long-term behaviors, with the initial conditions determining which long-term behavior is seen. This means that $\bar{Y}$ does not have a unique value for that parameter set.

4.   Simulations and sensitivity analysis

4.1.   Simulations

Some simulations for the monotone example function $\omega(Y) = \Sigma \left(1- e^{-aY} \right)$ are shown in Figure 9. The simulations use typical values for most parameters, as explained in Section 2.4, with $a = 2$ and two different values for the population flexibility parameter $\Sigma$ (see the open circles in Figure 5). The initial values of $Y$ were taken to be 50% and 20% higher than equilibrium for $\Sigma = 3$ and $\Sigma = 10$, respectively, while initial conditions for all other variables were chosen to be the relevant equilibrium values.

As predicted, the parameter set with smaller $\Sigma$ results in a stable EDE with $Y = 1.037$, while the larger $\Sigma$ results in an unstable EDE. We see from the simulation that the latter case produces oscillations.

Panel d shows additional detail that indicates how the oscillation occurs. Rapid changes in $Y$ cause rapid changes in $\omega$, which result in changes in $W$ that lag comparatively behind. At the peak in $Y$, $\omega$ is large and $W$ increases, which causes $Y$ to decrease. As $\omega$ subsequently decreases, along with $Y$, $W$ reaches its peak and then decreases in turn. Note that the amplitudes of oscillation are large for $Y$ and $\omega$, but small for $W$, owing to the differences in time scales on which these quantities change.

The horizontal black line in panel a shows the value of $Y^*$ for the oscillatory ($\Sigma = 10$) case. The peaks of the oscillation are far from the equilibrium value, which is why the mean value $\bar{Y}$ is significantly larger than $Y^*$, as previously seen in Figure 8.

For the non-monotone case, we consider a simulation designed to verify the features indicated in Figure 7. The dotted curves in panels a–c of Figure 10 show the case corresponding to the leftmost bullet point marked in Figure 7b, while the solid curves correspond to the rightmost bullet point ($d = 0.65$ and $d = 0.8$, respectively). From this plot, there should be three equilibria for each case, with the second unstable and the third stable. The chosen initial conditions are close to the small $Y$ equilibria, and the simulation plots show convergence to the equilibrium for $d = 0.65$ and the emergence of long-term oscillations for $d = 0.8$, as predicted by the theoretical results. Panel d shows behavior that matches that for the corresponding oscillatory case in Figure 9.

4.2.   Sensitivity analysis via the Morris method

The general trend seen in Figure 8 suggests that $\bar{Y}$ is less sensitive to $a$ and $\mathcal{R}_0$ than is $Y^*$, while both seem very sensitive to $\Sigma$ and less sensitive to $h$. However, this plot only looks at data in a small portion of the parameter space. To better understand how the mean infectious population is affected by the model parameters, we conducted a global sensitivity analysis using the Morris method ^[24]. We give a brief explanation of this method and an interpretation of the results for our model.

The Morris method takes a function $f(t; p_1, \dotsc, p_n)$ of one variable and $n$ parameters, and determines a distribution $F_i$ of elementary effects for each parameter, which we can assume are uniformly distributed on $(0, 1)$. One can use the inverse CDF to work with parameters not uniformly distributed on $(0, 1)$. To calculate $F_i$, first produce a sample $S$ from the parameter space $(0, 1)^n$. For each sampled $s \in S$, compute $s + e_i \Delta$ where $\Delta$ is a predetermined perturbation size and $e_i$ is the coordinate of the $i$th parameter. In other words, $s+e_i\Delta$ is obtained from $s$ by perturbing the $p_i$-value by $\Delta$. Also, for each $s \in S$, we compute the elementary effect $F_i (s) = (f(s + e_i \Delta) - f(s)) / \Delta$. This elementary effect, or simply effect, of $p_i$ at $s$ is a measure of $f$'s sensitivity to $p_i$. This process provides $n$ distributions $F_i$, each of size $|S|$. We can then compare the means and standard deviations of each distribution. If the mean of distribution $F_j$ is larger than that of $F_k$, then that would indicate that parameter $p_j$ has a larger effect on $f$ than $p_k$, on average across $S$. If the standard deviation of distribution $F_j$ is larger than that of $F_k$, then that would indicate that the effect of $p_j$ varies more throughout $S$ than the effect of $p_k$ does. However, one shortcoming of taking the mean of the elementary effect distributions is that of Type Ⅱ errors, namely, a mixture of perturbations of $p_i$, where some increase $f$ and some decrease $f$, can decrease the distribution mean, falsely indicating a small effect. To correct for this, we compute the distributions $G_i = |F_i|$ by simply considering the magnitude of changes in $f$.

We apply this method to the long-term mean prevalence $\bar{Y}$, with a parameter space built from the uniform distributions $\Sigma \sim U(2, 20)$, $R_0 \sim U(2, 10)$, $a \sim U(1, 5)$, and $h \sim U(0, 40)$, and with $v = 0.5$.³ From the parameter space, we took $10,000$ samples and perturbed each using $\Delta = 0.01$, to which we applied the inverse CDF prior to building the elementary effect distributions $F$ and $G$ for each parameter. The $F$ distributions for each parameter are displayed in Figure 11, and the statistics of the $F$ and $G$ distributions, rounded to two decimal places, are presented in Table 5.

³The value of $v$ does not affect the equilibrium solution, but when the equilibrium is unstable, it does change the infectious mean.

We see that $\mathcal{R}_0$ and $h$ have the largest and smallest effect on $\left| \bar{Y} \right|$ on average, respectively. The standard deviations of the $F$-distributions indicate that the effect of $\mathcal{R}_0$ varies the most across the parameter space, whereas the effect of $h$ varies the least. The distributions have long tails that skew the mean effect size, most dramatically for $\Sigma$, whose mean is clearly pulled to the left of its mode. Note that these long tails are not depicted in Figure 11, which truncates effect sizes whose magnitude is above 1. Finally, the discrepancy in magnitude between the $F$ and $G$ means indicates that $\bar{Y}$ is not monotonic in any parameter.

5.   Discussion and conclusions

In this study, we created and analyzed a relatively simple model that incorporates attitude change toward vaccination in response to the current disease prevalence. Asymptotic approximation allowed us to obtain simple conditions for the existence and stability of disease-free and endemic disease equilibria. We then chose some functional forms for the rate constants in the attitude change processes to provide examples, both for the more rational case of monotone increasing response with increasing prevalence and the case where a maximum response occurs at an intermediate prevalence level, with a subsequent decrease owing to vaccine fatigue.

5.1.   Effects of attitude change

Much of the work reported here focuses on instability in the model. The precursor model, which takes attitude to be immutable, does not show any instability ^[21]. The results suggest that attitude change is an important factor in the population dynamics of infectious diseases when vaccination is a part of the public health strategy.

5.2.   Comparison with the information index approach

The closest prior work to ours is the collection of papers that use an information index to modify the vaccination rate ^{[10,11,12,13,14,15,16]}. The two post-COVID papers in this group ^[15,16] lack a complete analysis of endemic disease equilibria. The earlier papers contain a complete EDE analysis; in each case, a larger sensitivity of vaccination to the information index can lead to an unstable EDE with a stable limit cycle. The primary differences between the information index models and our flexible attitude model are the mechanism by which changes in prevalence impact vaccination rates and the incorporation of loss of immunity in our model.

Even when the information index is based on recent prevalence history rather than just current prevalence, vaccination rates based on it respond much faster to prevalence changes than they do in our model. This is because of the different assumptions about human behavior. Prior to COVID-19, it was natural to assume that decisions about vaccination were based on a rational comparison of the payoffs of vaccination as compared to non-vaccination. Rational choices happen quickly, on the time scale of the disease processes. Following COVID-19, it is clear that vaccination decisions are based more on individual attitudes toward vaccination rather than rational comparison of alternative strategies. Unlike rational assessment, attitude changes are slow because people tend to be resistant to changing their opinions in response to new facts. In our model, vaccination rates respond to changes in prevalence on a slower time scale because the mechanism is indirect, with prevalence determining rates of attitude change rather than attitude itself and rate changes occurring on a slower time scale than that of prevalence changes.

Another implication of the difference between rational choice and attitude is that it makes sense to prescribe a mechanistic model for rational choice, whereas mechanistic models are less clear for human behavior. While it would be natural to expect an increase in a pro-vaccination attitude as disease prevalence increases, we also considered the possibility of a non-monotone response, in which a sufficiently high prevalence eventually results in a loss of confidence in the vaccine, thereby resulting in a decrease in pro-vaccination attitude with increasing prevalence.

For the monotone case, there is always a unique endemic disease equilibrium whenever the disease-free equilibrium is unstable. This unique EDE is sometimes stable, but when the response of the rate change coefficient to disease prevalence is strong enough, the EDE becomes unstable and the resulting behavior converges to a limit cycle.

The non-monotone case always has at least one EDE, but it can have up to three in the example function family we used, and possibly more for an unrealistically complicated function family. As seen in Figure 7, a variety of possibilities for stability exist, including either a stable EDE or a stable limit cycle when there is a unique EDE, and, when there are three EDE's, either a stable small-$Y$ EDE or a stable limit cycle, in addition to an unstable medium-$Y$ EDE and a stable large-$Y$ EDE.

Of particular interest are the amplitudes of the Willing subpopulation and the Unwilling to Willing transition rate, as shown in the plots of $\omega/\Sigma$ in Figures 9b and 10b. The equilibrium for the stable example in Figure 9 is 66% Willing (dotted curve in panel c), but the dotted curve in panel b shows that attainment of this level of willingness requires that the Unwilling to Willing transition coefficient ($\omega$) be much larger than the reverse transition coefficient ($\Sigma-\omega$). This is because newborns arrive in the population as Unwilling and become Willing only through the transition.

5.3.   The $\epsilon \to 0$ approximation

The use of asymptotic methods to produce approximate analytical results has been studied in detail by Ledder ^[22], with the conclusion that the actual differences in prediction between results for $\epsilon \to 0$ and small values like $\epsilon = 5\times 10^{-4}$ are often insignificant compared to differences caused by small changes in the estimates used for the model parameters. This claim can be tested for examples in this study.

Using the asymptotic limit $\epsilon \to 0$ simplified the mathematical analysis by reducing the system from five variables to four. We could instead have done the eigenvalue analysis for the full system, which would yield stability boundaries without having simple analytical results. The details appear in Appendix C.

Figure 12 shows the same scenario as in the right panel of Figure 5 but comparing the results for $\epsilon = 5\times 10^{-4}$ and those for $\epsilon \to 0$. As before, the heavy blue curve is the stability boundary calculated with the $\epsilon \to 0$ approximation. For comparison, the stability boundary for the 5-component system is the dashed black curve. The difference is noticeable, but small. If trying to fit actual data, it would be worthwhile to use the more accurate result; however, our work uses a generic model rather than one for a specific disease. Also, the curves are far more dependent on the other parameter values than on the difference between $\epsilon = 5\times 10^{-4}$ and $\epsilon = 0$.

5.4.   Conclusions

While the model is speculative in nature, it does suggest the possibility that attitude changes in response to prevalence can destabilize an endemic disease equilibrium, with multiple long-term outcomes in the case where the response is non-monotone. It would perhaps be difficult to identify this type of instability in real data because it would have seasonal and stochastic variation superimposed on it, but it would be worth exploring.

In searching for real examples of EDE destabilization, it is helpful to know what disease features are likely to promote destabilization. Most obviously, the sensitivity of the rate constants for attitude changes to differences in prevalence should be key, as evidenced by the greater tendency toward instability in the parameters $\Sigma$, $a$, and $d$ in the functions we used for these rate constants. It is also worth noting the destabilizing influence of the rate constant for loss of immunity $h$. While we might expect a faster loss of immunity to be destabilizing, the plots of Figure 5 suggest otherwise, as they show that a larger value of the attitude variability parameter $\Sigma$ is needed for instability for a relatively large value of $h$ as compared to a disease where there is no loss of immunity. But note also from this figure that the overall incidence of the disease is much higher for diseases with short-lived immunity than diseases with more permanent immunity, as of course would be expected.

While it is customary to consider the equilibrium values as measures of disease impact, our work shows that this is only true when the endemic disease equilibrium is stable. When there is a stable limit cycle, the mean prevalence is more important and can differ noticeably from the equilibrium prevalence.

In an effort to focus on the process of attitude change motivated strictly by disease incidence, we made a number of simplifying assumptions: simple SIR disease class structure, perfect vaccination that prevents infection and (through boosters) loss of immunity, no influence of incidence history, and no drivers of attitude connected to interpersonal or media-based communication. Of these assumptions, the assumption of perfect vaccination is the most likely to make a qualitative difference in outcomes, as imperfect vaccination would make prevalence less sensitive to vaccination rates and likely reduce the tendency toward instability.

While our study does not directly apply to COVID-19, it does suggest possibilities for how attitude changes will impact COVID-19 prevalence. While the actual epidemiological structure for COVID-19 should include a latent (usually called 'exposed') class and an asymptomatic class, these differences are unlikely to affect qualitative outcomes in the long run. Gradual loss of immunity is particularly important for COVID-19 population dynamics, with a typical loss of immunity over 6 months corresponding to a much larger value of $h$ than we used in our generic disease examples.

Use of AI tools declaration

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

Acknowledgments

This material is based on work done at the American Institute of Mathematics (AIM) with support from the National Science Foundation (NSF) under grant 2015462. The authors thank both institutions for their support. Lee gratefully acknowledges the financial support of Stonehill College (2023 Professional Development Grant).

Conflict of interest

Glenn Ledder is a guest editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. The authors declare there is no conflict of interest.

Appendix

A.   Proposition 1

Given $Y = 0$ for the disease-free equilibrium, the $W$ and $P$ equations have equilibrium relations

where the parameters $\psi$ and $\omega$ and those that depend on them are all understood to be evaluated at $Y = 0$. Substituting these relations for $W$ and $P$ into the equilibrium $S$ equation yields the disease-free equilibrium susceptible fraction

The vaccine-reduced reproduction number $\mathcal{R}_v$ can be found from the scaled problem (2.8) using first principles. If a small population $Y_0$ of infectious persons is introduced into a population that is at the disease-free equilibrium, the rate of infections in the scaled model from that initial attack will be $\mathcal{R}_0 S^* Y_0$, and the rate of removal of those infectious individuals will be $Y_0$. The average number of new infections per initial infectious individual will be the ratio of these two rates, that is

To show the connection between $\mathcal{R}_v$ and stability, we begin with the Jacobian, evaluated at the DFE:

where $\Gamma = \epsilon^{-1}$ and all state variables are evaluated at the equilibrium and all functions of $Y$ at $Y = 0$. Lower case letters are defined in (2.5). The block structure of the matrix allows us to immediately identify two of the eigenvalues as

The remaining eigenvalues come from the center block. Since this block has negative trace and positive determinant, the associated eigenvalues have negative real part; hence, stability is determined entirely by $\mathcal{R}_v < 1$.

From (A.2), rearranging $\mathcal{R}_v < 1$ yields a minimum requirement for $\omega$, given all other parameter values,

where $\bar{\Sigma}$ is evaluated at $Y = 0$ and $r = \mathcal{R}_0-1$.

B.   Proposition 3

The Jacobian matrix at the endemic disease equilibrium is

where

have been introduced to simplify the notation.

The usual way to compute the characteristic polynomial for a matrix $J$ is to find the determinant of $\lambda I-J$; however, it is more efficient to use the characteristic polynomial theorem ^[25]:

Theorem 1. For an $n \times n$ matrix $J$, let $S$ be the set of all nonempty subsets of the integers $1, 2, \ldots, n$. For each possible $K \in S$, let $J_K$ be the determinant of the principal submatrix of $J$ that contains the entries in the rows and columns indicated by the index set $K$. Then the characteristic polynomial of $J$ is

where

For (B.1), we can quickly get leading order results⁴

⁴Because $\epsilon \to 0$, we have $\Gamma \to \infty$; thus, the term in each $c_m$ with the most factors of $\Gamma$ dominates the other terms in the formula for that $c_m$.

To leading order, we then have the characteristic polynomial as⁵

⁵The reader may wonder why we must keep all the terms in this polynomial even though some have more factors of $\Gamma$ than others. The answer is that the magnitude of $\lambda$ is different for different roots. For example, the assumption $\lambda = O(1)$ yields a linear polynomial to leading order. This correctly approximates only one of the four eigenvalues; the other three must have $\lambda \to \infty$ as $\epsilon \to 0$.

where

Using only the leading order terms, we obtain the Routh array ^[22] as

where

Stability requires the quantities in the first column of the array to be positive. For simplicity, we can substitute $k_3 > 0$ for the more complicated $q_1 > 0$. To see this, suppose $k_4 > 0$ and $q_2 > 0$ as required. The formula for $q_2$ then implies that $k_3$ and $q_1$ must have the same sign. Stability requires $q_1$ to be positive, so it requires $k_3$ to be positive also. But then $k_3 > 0$, along with $k_4 > 0$ and $q_2 > 0$ is sufficient to guarantee the required condition $q_1 > 0$.

After simplification, the three requirements reduce to

Taking $\Sigma' = 0$ simplifies the stability criteria. Suppose $\omega' > 0$, as is the case for monotone functions and for $Y$ small enough for non-monotone functions. Then $A > 0$, so the first two conditions are automatically satisfied. If $\omega' < 0$, then the limiting criterion is always $k_4 > 0$, or $A > - \bar{\Sigma}/\rho = {\omega^*}'$. If this condition is satisfied, then the first condition is satisfied as well; hence, $1+A > 0$. The third condition is now satisfied automatically because both factors on the left side of the inequality are positive while the right side of the inequality is $\mathcal{R}_0 hvA < 0$. Taking both cases into account, the stability criteria simplify to

where

C.   Details for the original five-component system

Here we provide the details needed to do numerical stability analysis for the endemic-disease equilibrium for the original five-component system (2.7).

From the $W$ and $P$ equations, we obtain the results

where $y = \mathcal{R}_0 Y$ and $p = \mathcal{R}_0 P$ as for the four-component system. The $Z$ equation then yields

whereupon the $S$ equation yields the result

Combining this result with $\omega = \xi p$, we see that solutions occur when

This generalizes the formula obtained for $\omega^*$ for the four-component case.

With the additional assumption $\Sigma' = 0$, we obtain the Jacobian

Now suppose we want to find the stability boundary to compare with the solid blue curves in Figures 5 and 7(c). Using the $Y$ value at equilibrium as a parameter, we first calculate the corresponding values of $p$ and $\omega^*$ from (C.2) and (C.3). Then we use $\omega = \omega^*$ to calculate the remaining parameter for $\omega(Y)$ [$a$ for $\omega_1$ and $d$ for $\omega_2$], and then use the analytical derivative formula for $\omega'$. We can then think of the Jacobian as being an explicit function of $Y$ and identify the stability boundary for given $\Sigma$ [and $a$ for the case of $\omega_2$] as the value for which the dominant eigenvalue is 0.