Stochastic games of parental vaccination decision making and bounded rationality

1.   Introduction

In disease modeling, applying game theory to understand human behavior about mitigation choices offers crucial insights into the dynamics of epidemics. Conventional methods often presuppose rational decision-making; however, as the COVID-19 pandemic has demonstrated, human behavior can sometimes hinder public health initiatives to control the spread of the disease. Individual reactions to pharmaceutical interventions such as vaccination and non-pharmaceutical measures, including mask use, are shaped by a complex mixture of factors that include social norms, misinformation, and trust in authorities ^[1].

Moreover, integrating concepts from behavioral economics and psychology, such as bounded rationality and prospect theory, adds realism to these models. These theories recognize that individuals do not always process information flawlessly or make decisions based solely on maximum utility ^[2]. Instead, decisions are influenced by limited cognitive resources, biases, and how information is framed.

Kermack and McKendrick (1927) initially laid the groundwork for deterministic models of epidemics, which have since evolved to include stochastic components ^[3]. In ^[4], a pioneer in stochastic disease modeling was introduced using stochastic processes like birth-and-death processes. In that study, the authors introduced stochastic elements into mathematical epidemiology, making it possible to model the inherent randomness of disease transmission and study extinction thresholds to better understand and predict the dynamics of infectious diseases. Several researchers have studied stochastic models of the spread of diseases in the form of stochastic differential Eqs (SDEs). SDEs incorporate randomness directly into the model, representing various sources of uncertainty and variability in disease transmission and progression. These models have been particularly beneficial in capturing the random nature of contact patterns and the influence of random events on the spread of infectious diseases. Incorporating SDEs into epidemic modeling provides a more comprehensive understanding of how diseases spread and how they can be controlled under uncertainty.

Allen (1994) provided foundational work on stochastic epidemic models, demonstrating how SDEs can model demographic stochasticity and environmental variability in disease dynamics ^[5]. Gray et al. (2011) extended this approach to the SIS model, exploring the existence of a stationary distribution and derived expressions for its mean and variance ^[6]. Tornatore et al. (2005) investigated the stability of a stochastic SIR system, providing insights into the conditions under which disease extinction or persistence occurs ^[7]. Furthermore, Gao et al. (2019) studied a delay differential Eq SIS model with a general non-linear incidence rate, highlighting the impact of delays and non-linearities in disease transmission ^[8]. Cao et al. (2017) focused on the dynamics of an epidemic model with vaccination in a noisy environment, emphasizing the role of stochastic effects in vaccination strategies ^[9]. Liu (2023) modeled multiple transmission pathways in a scenario of waterborne pathogens using a stochastic model, illustrating the complexity and variability in disease transmission routes ^[10]. See also ^[11] for an introduction to stochastic epidemic models.

Researchers have continued to advance the field by exploring the effects of multiple noise sources on infectious disease models. For example, Jian et al. (2024) examined the impact of multisource noise on disease dynamics, offering insights into how different environmental randomness can influence disease spread ^[12]. Babaei et al. (2023) applied SDEs to model the COVID-19 pandemic, integrating stochastic elements into SEIAQHR models to capture the variability and uncertainty of real-world scenarios ^[13]. Similarly, Iddrisu et al. (2023) highlighted the critical role of stochasticity in understanding cholera epidemiology, demonstrating how SDEs can provide more accurate predictions and inform public health interventions ^[14]. These contributions underscore the importance of stochastic models in epidemiology, as they provide more realistic and flexible frameworks for predicting disease outcomes and formulating effective public health policies.

Rational decision-making theories often disregard information processing constraints, presuming that rational decision-makers always opt for the highest rewards without factoring in the effort or resources needed to determine the best action ^[15,16,17]. H.A. Simon's bounded rationality theory addresses these deviations, suggesting that decision-makers operate under limited information processing resources ^[18,19]. Behavioral economics experiments reveal deviations between perfectly rational decision-makers and those constrained by bounded rationality ^[20].

Behavioral epidemiology has evolved in the last two decades to become a major domain in the research of infectious diseases. The intricate interplay between human behavior and epidemiological outcomes highlights this. Foundational contributions to the literature of behavioral epidemiology, such as by ^{[21,22,23,24,25]}, have shaped the understanding of the impact of behavioral responses to disease outbreaks and their control. In particular, it gave insight into how parental vaccination decisions affect childhood disease dynamics and public health interventions, especially with vaccine scare ^[26,27,28]. Free-riding behaviors ^[29], social networks ^[30], social learning ^[31], and economic incentives ^[32], to name a few factors, have been shown to impact vaccination uptake and thus disease spread. Those research findings emphasize the importance of integrating behavioral and epidemiological modeling approaches to address challenges to disease control ^[25]. They demonstrate the complex interplay of behavioral, social, cognitive, and economic factors in human behavior in general and vaccine decision-making, in particular, which form the basis for this research.

We seek to add a building block to behavioral epidemiology advances regarding factors influencing parental vaccination decisions, particularly rationality in the decision-making process. We extend the work in ^[28], where the authors emphasized integrating social factors into behavioral responses and infectious disease models, and ^[27], where the authors explored the use of prospect theory to enhance replicator dynamics in vaccination decision making. Here, we aim to build on these models by focusing on bounded rationality in decision-making resulting from stochasticity in the perception of the payoff of vaccination, leading to a mutation term in the replicator dynamics. Our results contribute to understanding parental behavior towards childhood immunization by showing how collective deviations in those perceptions can make disease persist. We also found that the consequences of these deviations could be counteracted using an optimal temporal discount in the cost of vaccination.

In this paper, we derive a stochastic differential Eq of RD for vaccination ^[33]. The stochasticity arises from the perception of the utilities of the choices and deviations from rationality. Using stochastic game theory to model vaccine choice represents a significant step forward in capturing this complexity. In this framework, a collective decision to vaccinate is influenced not only by personal perceptions of risk-benefit but also by the choices made by others in their community and external but common environmental, cultural, informational, and economic factors. These shared factors are essential for understanding collective dynamics and population-level responses under systemic uncertainty. This interdependence of decisions, along with the inherent variability in the way individuals assess risks and benefits under uncertainty, can also be effectively modeled as a stochastic game in which the outcomes depend randomly on the actions of all participants ^[28,34,35]. To our knowledge, we are the first to introduce and study a stochastic behavioral game theory model in the disease modeling literature.

The paper is summarized as follows. In Section 2, we derive a new stochastic behavioral game theory Eq that models vaccination versus non-vaccination choices. We pair that Eq with a stochastic SIR model in which the transmission rate also evolves randomly. This will be followed in Section 3 by an analytical stability analysis of the model's equilibria. We also use stochastic simulations of stochastic differential Eqs to further study the stability of the model equilibria and corroborate some of the analytical results. Additionally, we use stochastic optimal control of the cost of vaccination to counteract the effect of bounded rationality. We follow that section with a discussion and conclusion section.

2.   Materials and methods

Let $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq 0}, \mathbb{P})$ be a complete probability space, with filtration $\{\mathcal{F}_t\}_{t\geq 0}$, with respect to which all Brownian motions $B_i$ and $W_i$ introduced in the following are defined. We use a stochastic replicator Eq following the methods of establishing a stochastic replicator Eq in ^[36,37] to guarantee the feasibility of solutions of the stochastic RD Eq. The stochastic RD Eq models the evolution over time of strategies or traits in a population by incorporating random fluctuations due to environmental variability or when the population is finite. For instance, let a population have two strategies or traits, A and B, like vaccination and non-vaccination strategies. If strategy A has a higher payoff than B, then the RD Eq's deterministic component would imply that the adopters of strategy A must increase. However, the stochastic component could cause temporary increases in the adopters of strategy B, while strategy A is more gainful and vice versa.

Let $X(t) = (X_1(t), X_2(t))$ be the number of vaccinators and non-vaccinators in a community of a total size of $N: = X_1+X_2$ decision makers and let the fractions be denoted by $x_i = X_i/N$. Since $x_1+x_2 = 1$, we will also use the lower case $x$ to denote $x_1$ whenever it does not cause confusion. Let $u_i(x_1, x_2) = \sum_{j = 1}^2 v(s_i, s_j)x_j$ be the expected payoff for a decision maker adopting strategy $s_i$ if randomly and uniformly matched with a decision maker adopting strategy $s_j$. The function $v(s_i, s_j)$ is the payoff for a decision maker who adopts strategy $s_i$ interacting with a decision maker adopting strategy $s_j$.

A deviation from rationality due to external and internal factors is modeled here by adding a white noise term to the payoff $v(s_i, s_j)$ giving a new payoff function $v(s_i, s_j)+\nu_i \dot{B}_i(t)$. The stochastic process $\dot{B}(t) = (\dot{B}_1(t), \dot{B}_2(t))$ is a two-dimensional white noise. Volatility $\nu_i > 0$ is the magnitude of the noise in the perceived utility. We assume here that all vaccinators and all non-vaccinators share the same volatility $\nu_1$ and $\nu_2$, respectively. Maintaining collective population dynamics, the stochastic RD model postulates that stochastic noise affects payoff perceptions uniformly across adopters of the same strategy. This does not account for individual variability in the perceived payoffs. That assumption reflects the influence of common external factors on risk perceptions, aligning with psychological theories of collective decision-making under uncertainty.

The following stochastic differential Eq then gives the replicator Eq modeling behavioral dynamics,

or

where $\kappa > 0$ is the replication rate. Applying Itô's formula (see Lemma S2 in the Supplementary) to the transformation $x_1 = X_1/N$, we get

Since $\nu_1 dB_1(t)-\nu_2 dB_2(t)$ has the same distribution as $\sqrt{\nu_1^2+\nu_2^2}\, dW(t)$, where $dW(t)$ is also Gaussian with mean zero and variance $dt$, then Eq (2.3) is equivalent to:

Eq (2.4) guarantees that $0\leq x_1(t)\leq 1$ for all $t > 0$ if $0\leq x_1(0)\leq 1$.

We model the payoffs of vaccination to be $v(s_1, s_2) = -\omega+\delta\, x_1$ and of vaccination to be $v(s_2, s_1) = -I+\delta\, x_2$, see ^[28]. In terms of the vaccination choice model, $\kappa$ is the social learning rate, $\omega$ is the cost of vaccination, and $\delta$ is the group pressure. Social norms of vaccination and non-vaccination enforced by vaccinators and non-vaccinators add $\delta\, x_1$ and $\delta\, x_2$ to the payoffs of vaccination and non-vaccination strategies, respectively.

The noise in perceived payoff results in the introduction of a new mutation term $\kappa^2 x_1 x_2 (\nu_2^2 x_2- \nu_1^2 x_1)$ to replicator dynamics; see, e.g., ^[36,38]. A mutation term induces irrational choices of a strategy that might not be optimal. In terms of vaccination, for certain large values of $\nu_1$, the mutation term is a bounded rational term that can cause agents to adopt the non-vaccination strategy while it has a smaller payoff; that is $v(s_2, s_1) < v(s_1, s_2)$.

The transmission dynamics are assumed to follow a susceptible-infected-recovered (SIR) compartmental model. Putting $x_1 = x$ and $x_2 = 1-x$ in Eq (2.4) and coupling it with the SIR model yields the system of Eqs:

where $S$, $I$, and $R$ are the proportion of susceptible, infected, and recovered individuals in the population at time $t$, such that $S+I+R = 1$. The parameter $\mu$ is the per capita birth/death rate, $\beta$ is the transmission rate, $\gamma$ is the recovery rate. Recall that $\kappa$ is the social learning rate, $\omega$ is the cost of vaccination, and $\delta$ is the group pressure. The term $\mu(1-x)$ is the recruitment rate of non-vaccinated children to the susceptible compartment due to parental choice.

Let $\sigma_1$ be the magnitude of noise in disease transmission. The last Eq in (2.5), a stochastic replicator Eq, follows the replicator model introduced by ^[37], where it is shown to have solutions $x$ between zero and one; see above. Notice that we changed the notation $\nu_1$ to $\sigma_2$ and $\nu_2$ to $\sigma_3$, to be in order within the coupled model. If $\sigma_2^2 = \sigma_3^2 = \frac{\delta}{\kappa}$ then, $\delta\, (2\, x-1)+\kappa \left(\sigma_3^2 - (\sigma_2^2+\sigma_3^2) x \right) = 0$. In this case, the collective noise would cancel out the group pressure in the stochastic replicator Eq. Hence, it seems that bounded rationality could nullify or counteract social norms.

Since $S+I+R = 1$, then the third Eq of $R$ in Eq (2.5) is redundant, so henceforth we will consider the model (2.5) without it and stop using the recovered compartment. Define the set of solutions for the model in Eq (2.5) as $\mathcal{S} = \{(s, i, x)\in \mathbb{R}_+^3: 0\leq s+i \leq 1 \text{ and } 0\leq x \leq 1\}$. In the following subsections, let $(S(t), I(t), x(t))$ be a solution of the model (2.5) with initial state $(S(0), I(0), x(0))\in \mathcal{S}$. The proofs of existence, uniqueness, and boundedness of the solution of the model (2.5) are given in S1.2.

3.   Results

In this section, we discuss the different equilibria of the model in (2.5) and their stability. We note that some specific states of the three stochastic processes $S$, $I$, and $x$ are absorbing states. The states $0$ and $1$ are absorbing states for $x$ and $0$ is an absorbing state for $I$.

We will support the results in this section with simulations using the Milstein algorithm for the numerical simulation of SDEs ^[39] implemented using Python. We will use parameter values from ^[28], such as $\mu = 1/50$ year$^{-1}$, $\gamma = 365/22$ year$^{-1}$, and $\kappa = 1.69$ year$^{-1}$. We have to mention here that most of the simulated processes exhibit slow mixing rates, and hence they take a very long time to converge to equilibrium.

3.1.   Model's equilibria

The model has five equilibria of the drift term (the deterministic part) in Eq (2.5). Those five equilibria and their existence conditions are found to be those of ^[28] when $\sigma_1^2 = \sigma_2^2 = \sigma_3^2 = 0$. Three of those are disease-free equilibria:

1. full vaccine uptake and no susceptibility, $\mathcal{E}_1 \equiv (S_1, I_1, x_1) = (0, 0, 1) \, ,$

2. no vaccine uptake and full susceptibility, $\mathcal{E}_2 = (1, 0, 0)$, and

3. partial vaccine uptake and partial susceptibility, $\mathcal{E}_3 \equiv (1-x_3, 0, x_3)$, where

The equilibrium $\mathcal{E}_3$ exists in two regions:

and

The two boundaries of either region intersect at $\left(\frac{\kappa}{2}(\sigma_2^2+\sigma_3^2), \frac{\kappa}{2}(\sigma_3^2-\sigma_2^2) \right)$.

The other two are endemic equilibria that exist when the basic reproduction number $R_0 = \dfrac{\beta}{\mu+\gamma} > 1$.

They are:

1. no vaccine coverage, $\mathcal{E}_4 \equiv \left(\frac{1}{R_0}, \frac{\mu}{\mu+\gamma}\left(1-\frac{1}{R_0}\right), 0\right)$, and

2. partial vaccine coverage, $\mathcal{E}_5 \equiv \left(\frac{1}{R_0}, \frac{\mu}{\mu+\gamma}\left(1-\frac{1}{R_0}-x_5\right), x_5\right),$ where

The equilibrium $\mathcal{E}_5$ exists in two regions

and

The two boundaries of either region intersect at $\left(\dfrac12 \dfrac{\mu}{\mu+\gamma}+\frac{\kappa}{2}(\sigma_2^2+\sigma_3^2), \dfrac12 \dfrac{\mu}{\mu+\gamma} (1-\frac{2}{R_0})+\frac{\kappa}{2}(\sigma_3^2-\sigma_2^2) \right)$.

The disease-endemic equilibria appear in the next stability analyses as types of behavior and not as exact values. The equilibria $\mathcal{E}_3$ and $\mathcal{E}_5$ represent the partial acceptance of the vaccine. However, in the latter, partial acceptance of the vaccine does not lead to the eradication of the disease.

3.1.1.   Dynamical regimes of stability

Define the mean-value function $\overline{Y}(T)$ of any function $Y(t)$ over $[0, T]$ to be

and let its limit be defined as

Let $Y_* = \liminf\limits_{t\to \infty} Y(t)$ and $Y^* = \limsup\limits_{t\to \infty} Y(t)$.

Next, we define modes of stability, such as exponential stability, and introduce a new mode of stability that we call logistic stability. Logistic stability helps in studying the impermanence of dynamic processes that take value in the unit interval $[0, 1]$. To our knowledge, this concept has not been used before in dynamical system analyses.

Definition 1 (Almost Sure Exponential Stability). We say that a process $Y(t)\in [0, \infty)$ has an almost sure exponentially stable equilibrium at $y = 0$ if for some $c > 0$,

almost surely.

Definition 2 (Almost Sure Logistic Stability). We say that a process $Y(t)\in [0, 1]$ has an almost sure logistically stable equilibrium at $y = 0$ and $1$ if for some $c > 0$,

almost surely.

The unit interval represents a natural domain for many real-life phenomena, such as proportions, probabilities, or normalized populations constrained to the unit interval $[0, 1]$. The concept of logistic stability serves the need to characterize the transient dynamics of those types of processes as well as their behavior in convergence towards stability. Other convergence behaviors like exponential or asymptotic stability provide some types of convergence to equilibria. Those types may not capture other systems where the dynamics of the processes are bound and/or demonstrate logistic-like growth or decay which tend to be slower than the exponential ones. In infectious disease modeling, logistic stability could capture the impermanence of disease states within the unit interval, like prevalence, when epidemics go through transient phenomena. Logistic stability could distinguish those transient dynamics from other forms, emphasizing steady-state behavior. The limiting behavior of $\frac{1}{t} \log \frac{Y(t)}{1-Y(t)}$ can quantify the rate at which a disease dies out or persists, providing a new metric for control measures.

Define the stochastic basic reproduction number

which is always smaller than or equal to $R_0$.

First, we will study the global stability of a disease-free equilibrium. There are several cases in which the disease-free equilibrium is stable. The first case depends on the stochastic basic reproduction number and the magnitude of the noise in transmission.

Theorem 1. If $R_0^s < 1$ and $\sigma_1^2\leq \beta$ or if $\frac{\sigma_1^2}{\beta} > \max\left(1, \dfrac{R_0}{2}\right)$, then a disease-free equilibrium is exponentially globally stable almost surely.

Proof. We will prove the theorem in two parts, by showing that the following two inequalities hold:

(C.Ⅰ) If $\frac{\sigma_1^2}{\beta} > \max\left(1, \dfrac{R_0}{2}\right)$, then

(C.Ⅱ) If $R_0^s < 1$ and $\sigma_1^2 \leq \beta$, then

Let $V(t, I) = \log(I)$, then by Itô's formula (see Lemma S2 in the Supplementary),

which by integration over the interval $(0, t)$ gives

where $M(t) = \sigma_1 \int_0^t S(u) dW_1(u)$ is a local continuous martingale with $M(0) = 0$. Since

then by the Strong Law of Large Numbers, $\limsup\limits_{t\to \infty} \dfrac{M(t)}{t} = 0$ almost surely, , see ^[40] and Lemma S3 in the Supplementary.

Now to prove (C.Ⅰ), when $\frac{\sigma_1^2}{\beta} > \max\left(1, \dfrac{R_0}{2}\right)$, note that the quadratic function

attains its maximum value of $\frac{\beta^2}{2\sigma_1^2}$ when $S = \frac{\beta}{\sigma_1^2}$, which is in the interval $(0, 1)$ if $\sigma_1^2 > \beta$, in which case

and so

From this (C.Ⅰ) follows as

To prove (C.Ⅱ), when $\sigma_1^2 \leq \beta$, consider the function

which attains a maximum value on the interval $[0, 1]$ at $S = 1$ with maximum value

Thus,

and so

From this (C.Ⅱ) follows as

Theorem 1 seems to be a classic conclusion; see, e.g., ^[41]; but it is shown here to hold even with the presence of a stochastic process $x(t)$ in the first Eq. It shows that disease extinction could still be achieved here even when $R_0 > 1$ as long as $R_0^s < 1$, see Figure 1(a). Another case in which the disease could be eradicated while $R_0 > 1$ is when $\frac{\sigma_1^2}{\beta} > \max\left(1, \dfrac{R_0}{2}\right)$ as shown in Figure 1(b). Even when $R_0^s > 1$; or specifically, when the conditions of Theorem 1 are not satisfied, eradication of the disease is possible through vaccination, as we will see next.

In ^[28], using a deterministic model, it was shown that while $R_0 > 1$, the disease could be eradicated as the uptake of the vaccine increases to the level of $100\%$ when the peer pressure of the vaccination group is greater than the cost of the vaccine ($\delta > \omega$). However, it was also shown that if the cost of vaccination and the peer pressure of the non-vaccinating group exceed the risk of the disease ($\omega+\delta > \frac{\mu}{\mu+\gamma}(1-\frac{1}{R_0})$), parents might skip vaccination and the disease becomes endemic; see also ^[27]. Moreover, both equilibria were shown to be bistable. Here, bounded rationality changes those findings in a way that noise can impede the roles played by the cost of vaccination and peer pressure in ^[28]. The following theorem shows how noises in the perceived utilities can change the vaccination selection, and consequently the disease dynamics.

Theorem 2. If $\bar{I}(t)\xrightarrow[]{t\to\infty} I_0$, then $x(t)$ has an almost sure logistic stable equilibrium at $0$ and $1$ under the following conditions:

1. If $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} > \max(I_0-\delta, \delta)$, then $x(t)\xrightarrow[]{t\to\infty} 0$ logistically almost surely, and

2. If $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} < \min(I_0-\delta, \delta)$, then $x(t)\xrightarrow[]{t\to\infty} 1$ logistically almost surely, and also, $I(t)\xrightarrow[]{t\to\infty} 0$ and $S(t)\xrightarrow[]{t\to\infty} 0$ $a.s.$ whether $R_0^s$ is more or less than one.

Proof. Using Itô's formula (see Lemma S2 in the Supplementary),

and similarly

Therefore,

Since $\frac{W_2(t)}{t}\xrightarrow[]{t\to\infty} 0$ almost surely due to the Strong Law of Large Numbers of Brownian motion, then

almost surely, where

If $x_0 = 0$, which is equivalent to $x(t)$ approaches $0$ almost surely, then $L(x_0, I_0) = -\kappa \frac{\sigma_2^2-\sigma_3^2}{2}-\omega-\delta+I_0$. If $x_0 = 1$ and $I_0 = 0$, which is equivalent to $x(t)$ approaches $1$ almost surely, then $L(x_0, I_0) = -\kappa \frac{\sigma_2^2-\sigma_3^2}{2}-\omega+\delta$. The first and second values of $L(x_0, I_0)$ could have opposite signs for the same parameter values as when $I_0-\delta < \kappa \frac{\sigma_2^2-\sigma_3^2}{2}+\omega < \delta$. Therefore, $x(t)$ approaches $0$ logistically almost surely if $-\kappa \frac{\sigma_2^2-\sigma_3^2}{2}-\omega-\delta+I_0 < 0$ and $-\kappa \frac{\sigma_2^2-\sigma_3^2}{2}-\omega+\delta < 0$. Meanwhile, $x(t)$ approaches $1$ logistically almost surely if $-\kappa \frac{\sigma_2^2-\sigma_3^2}{2}-\omega+\delta > 0$ and $-\kappa \frac{\sigma_2^2-\sigma_3^2}{2}-\omega-\delta+I_0 > 0$.

An overlap region appears in Figure 2(a) as stated in the proof of Theorem 2 when $I_0-\delta < \kappa \frac{\sigma_2^2-\sigma_3^2}{2}+\omega < \delta$. In that region, both equilibria could be stable, as the simulation demonstrates below. The dashed region with vertical lines in Figure 2(a) is for $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} > I_0-\delta$ in which the disease-endemic and no-vaccination equilibrium $\mathcal{E}_4$-type equilibrium is stable; see Figure 2(b) for a simulation of a case in that region. Meanwhile, the dashed region with diagonal lines in Figure 2(a) is for $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} < \delta$, in which disease-free and full-vaccination equilibrium $\mathcal{E}_1$ is stable; see Figure 2(c) for a simulation of a case in that region. Again, in the common region, both $\mathcal{E}_1$ and $\mathcal{E}_4$ are stable, see below.

When $\sigma_2^2 = \sigma_3^2 = 0$, the two conditions of Theorem 2 are shown to be sufficient for the local asymptotic stability of $\mathcal{E}_4$ and $\mathcal{E}_1$, respectively, for the deterministic model in ^[28]. In particular, the latter equilibrium of the disease-free state was shown to be stable when $\delta > \omega$. However, the magnitudes of the white noise in the utilities of both groups, $\sigma_2^2$ and $\sigma_3^2$, influence the opinion and uptake of vaccination. In Figure 2(b), the disease cannot be eradicated even when $\delta > \omega$ is mainly due to the level of noise in the vaccinator group. Figure 2(c) represents a case where the conditions of Theorem 1 are not satisfied, but disease eradication is still possible as implied by part 2 of Theorem 2. If $R_0^s > 1$, the disease could be eradicated by increasing the peer pressure of the vaccination group or decreasing the noise in the utility of the vaccinator's group conducing to $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} < \delta$. Notice also that an increase in the noise of the non-vaccinator group could also increase vaccine uptake.

Achieving the goal of disease eradication in the vaccination campaign requires minimizing variability in decision-making and promoting consistent, rational choices throughout the population. However, that aim might not be plausible due to external factors. The noise in the perceived utility of the vaccination group could increase further, exacerbating disease eradication. When the social and stochastic external factors that influence utility functions are outside of the control of the policymaker, a policy to change the cost of vaccination is warranted. A stochastic optimal control approach (see e.g., ^[42]) could provide a policy for preventing infectious diseases by finding the appropriate cost of vaccination. As such, let the control variable $0 < u(t)\leq 1$ be the degree of reduction in the cost of vaccination. Including that function in the stochastic replicator Eq gives rise to the following Eq

The details of the stochastic optimal control are given in S1.3. Figure 3(a) provides a simulation in which $\mathcal{E}_4$ is stable. Full uptake of the vaccine and disease eradication is achieved, see Figure 3(b), with an optimal control $u^*(t)$ that reduces the cost of vaccination, see Figure 3(c). We note that the optimal reduction started to decline smoothly after being a constant over a relatively short period of time.

It would be difficult to control the disease without that discount in the cost of vaccination, especially with large noise in the vaccination group's utility function. The reason is in two scenarios: the first scenario is when $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} > I_0-\delta$ and $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} > \delta$, in which vaccination is completely rejected and the disease becomes endemic; see Figure 2(b). The second scenario is where $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} > I_0-\delta$ and $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} < \delta$, in which case both disease eradication and disease endemicity become stationary; see Figure 4. In Figure 4(a) the disease dies, while in Figure 4(b) the disease persists while all the values of the parameters and initial conditions are the same.

As shown in Figure 4, when both parts of Theorem 2 are true, namely $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} > I_0-\delta$ and $\omega+\kappa\frac{\sigma_2^2-\sigma_3^2}{2} < \delta$, then the bistability of both equilibria $x = 0$ and $x = 1$ (impermanence) occurs. A further examination of the stability of both equilibria was performed using a simulation study. Figure 5 shows the proportion of times $x(t)$ converges to $0$ and to $1$ in $200$ stochastic simulations of the model at different values of $\sigma_2^2$ and $\sigma_3^2$ with initial values $x(0) = 0.1, 0.2, \ldots, 0.9$. The bistability region shows a Bernoulli stationary distribution whose probabilities depend on the noises in the perceived utilities $\sigma_2^2$ and $\sigma_3^2$, as well as the initial acceptance of vaccination $x(0)$. The size of the bistability region could be seen to depend on $x(0)$ as shown through the panels of Figure 5. We notice that as initial acceptance of the vaccine $x(0)$ increases, the region in which full vaccination is stable in the $\sigma_2^2-\sigma_3^2$ plane expands. In other words, noises are more detrimental to vaccination campaigns if initial vaccine acceptance is low.

The following two theorems are generalizations to the disease dynamics in ^[28]. The following theorem identifies the boundaries of the limiting temporal means of processes when $R_0^s > 1$. It also shows the case when the disease cannot be eradicated at a low vaccination level; namely when its temporal mean of vaccine acceptance/uptake is below $1-\frac{1}{R_0^s}$.

Theorem 3. If $R_0^s > 1$ and $\bar{x}(t)\to x_0$, then a disease-endemic equilibrium will persist in the temporal mean. That is

1. For all $x_0\in [0, 1]$,

and

2. If $x_0 < 1-\frac{1}{R_0^s}$, then

and

Eq (3.13) holds true also when $x_0\geq 1-\frac{1}{R_0^s}$ from the upper inequality in Part (1).

3. If $x_0 < \frac{\beta}{\beta-\sigma_1^2}(1-\frac{1}{R_0^s})$ and $\beta > \sigma_1^2$ then

Proof. Since

then

where $\phi(t) = \frac{1}{\mu t}\left(S(t)+I(t)-S(0)-I(0) \right)$. Notice that $\lim_{t\to \infty} \phi(t) = 0$ $a.s.$ Since, $\bar{S}(t)\leq 1-\phi(t)-\bar{x}(t)$ for all $t > 0$. Therefore,

From Eq (2.5)

Then,

where $M_3(t) = \int_0^t \sigma_1 S(u)\, I(u)\, d{W}_1(u)$. Notice that $M_3(t)$ is a local continuous martingale with $M_3(0) = 0$ and $E(M_3(t)^2) = \sigma_1^2 \int_0^t (S(u) I(u))^2 du\leq \sigma_1^2 t$. Since $\limsup\nolimits_{t\to \infty} \dfrac{E(M_3(t)^2)}{t}\leq \sigma_1^2 < \infty$, then by the Strong Law of Large Numbers, $\limsup\nolimits_{t\to \infty} \dfrac{M_3(t)}{t} = 0$ almost surely, see ^[40] and Lemma S3 in the Supplementary. But since $x_0 < 1$, and $S(t)-S(0)\leq 1$ then (3.11) follows.

To prove part (2), note that

Thus

and so

Thus by Lemma S4 in the Supplementary,

if $x_0 < 1-\frac{1}{R_0^s}$. In which case, since

Then

Notice that $\frac{1}{R_0^s} > \frac{1}{R_0}$ so it is more than its deterministic counterpart.

To prove part (3), similarly, if $x_0 < \frac{\beta}{\beta-\sigma_1^2}(1-\frac{1}{R_0^s})$ and $\beta > \sigma_1^2$ then

Again,

and by Eq (3.15)

where $M(t) = \sigma_1 \int_0^t S(u) dW_1(u)$ is a local continuous martingale, with $M(0) = 0$, see ^[40] and Lemma S3 in the Supplementary. Since, $(\bar{S}(t))^2\leq \overline{S^2}(t)$, then

where

and

If we show that $\lim_{t\to \infty}\Xi(t) = \beta(1-\frac{1}{R_0^s})-(\beta-\sigma_1^2)x_0$ $a.s.$, then Lemma S4 in the Supplementary will imply that

if $x_0 < \frac{\beta}{\beta-\sigma_1^2}(1-\frac{1}{R_0^s})$. But $\lim_{t\to \infty}\frac{1}{t}\log(I(0)) = 0$ $a.s.$ and $\lim_{t\to \infty}\frac{1}{t}M(t) = 0$ $a.s.$ Finally,

and $\lim_{t\to \infty}\phi(t) = 0$ $a.s.$ complete this part.

In ^[28], where the deterministic version of the model is studied, the partial vaccination and persistence of the disease equilibrium $\mathcal{E}_5' \equiv \left(\frac{1}{R_0}, \frac{\mu}{\mu+\gamma}\left(1-\frac{1}{R_0}-x_5'\right), x_5'\right)$, where

is locally asymptotically stable when

The left hand side of inequality (3.16) is equivalent to $x_5' < 1-\frac{1}{R_0}$, the existence condition of $\mathcal{E}_5'$. Part 2 of Theorem 3, exhibits a similar result in which if the temporal mean of $x(t)$ converges to $x_0 < 1-\frac{1}{R_0^s}$, then the disease persists as in the deterministic model even when $x_0 > 0$. In which case, parts 2 and 3 of Theorem 3 give

Moreover, if $\sigma_1^2 = 0$ while $R_0 > 1$, then

When $x_0 = 0$, the disease persists in a way similar to $\mathcal{E}_4$.

When $R_0^s > 1$ some non-deterministic limits of the system occur under the condition that $\dfrac{\sigma_1^2}{\beta} < \dfrac{R_0}{2}$; see, e.g., ^[43]. Even when $\sigma_1^2 = 0$ and $R_0 > 1$, in view of Eq (3.18), disease prevalence still fluctuates; see Figure 6(a). The following theorem provides boundaries for those stochastic fluctuations.

Theorem 4. If $R_0^s > 1$, and $\dfrac{\sigma_1^2}{\beta} < \dfrac{R_0}{2}$, then

1. $\liminf\nolimits_{t\to \infty} S(t)\leq s_d \leq \limsup\nolimits_{t\to \infty} S(t)\leq 1-x_*$ a.s.

2. $\liminf\nolimits_{t\to \infty} I(t)\leq (1-s_d-x_*)\frac{\mu}{\mu+\gamma}$ a.s.

3. $(1-s_d-x^*)\frac{\mu}{\mu+\gamma} \leq \limsup\nolimits_{t\to \infty} I(t) \leq (1-s_d)$ a.s.

where $s_d: = \frac{1}{R_0}\, \frac{2}{1+\sqrt{1-\frac{2\sigma_1^2}{\beta R_0}}}\geq \frac{1}{R_0}$. Note that the first part implies that $x_*\leq 1- s_d$.

Proof. Since,

Let $\phi(z) = -\frac12 \sigma_1^2 z^2+\beta z-(\mu+\gamma)$, which has one zero $s_d$ between zero and one if and only if $\frac{\sigma_1^2}{\beta} < \dfrac{R_0}{2}$ and $R_0^s > 1$. The function $\phi(z)$ can have a peak at $z^* = \frac{\beta}{\sigma_1^2}$ such that $s_d < z^* \leq 1$ if and only if $\beta\leq \sigma_1^2$. Since $\phi(z)$ is increasing for $z < z^*$ then there exists $\epsilon > 0$ such that $\phi(z)$ is increasing for $z < s_d+\epsilon$. Thus, $\phi(z) < \phi(s_d-\epsilon) < 0$ for $z < s_d-\epsilon$.

Assume that $\limsup\nolimits_{t\to \infty} S(t) < s_d-\epsilon$ on a set $\Omega_1$ with positive probability, and so $0\leq S(u) < s_d-\epsilon < 1$ for all $u > t_0(\omega)$ for each $\omega \in \Omega_1$. Thus, $\phi(S(u)) < \phi(s_d-\epsilon) < 0$ for all $u > t_0(\omega)$ for each $\omega \in \Omega_1$. Then, $\limsup\nolimits_{t\to \infty} \frac{\log(I(t))}{t}\leq \phi(s_d-\epsilon) < 0$ with positive probability. That means $\lim_{t\to \infty} \bar{I}(t) = 0$ with positive probability, but that would contradict Theorem 3 part (2).

Now, assume that $\liminf\nolimits_{t\to \infty} S(t) > s_d+\epsilon$ on a set $\Omega_2$ with positive probability, and so $s_d+\epsilon\leq S(u) \leq 1$ for all $u > t_0(\omega)$ for each $\omega \in \Omega_2$. Thus, $\phi(S(u)) > \phi(s_d+\epsilon) > 0$ for all $u > t_0(\omega)$ for each $\omega \in \Omega_2$. Then, $\liminf\nolimits_{t\to \infty} \frac{\log(I(t))}{t}\geq \phi(s_d+\epsilon) > 0$ with positive probability. But that would contradict with that $\frac{\log(I(t))}{t} < 0$ for all $t$ since $I(t)\leq 1$. Thus, part (1) follows.

Again, since

then

Therefore, by the generalized L'Hôpital rule (see Supplementary)

and so

Also, by the generalized L'Hôpital rule (see Supplementary)

and so

if $x^* = \limsup\nolimits_{t\to \infty} x(t) < 1$. Since,

for all $t > 0$, then

Note that $s_d \downarrow \frac{1}{R_0}$ as $\sigma_1\to 0^+$. Figure 6(a) shows a simulation of one instance of the case where $\sigma_1^2 = 0$ and $R_0^s = R_0 = 6$, in which the proportion of infected fluctuates due to the stochastic nature of vaccine uptake. In that case, $S^* = S_* = \frac{1}{R_0}$. Figure 6(b) illustrates the case where $\sigma_1^2 > 0$ and $R_0^s = 4.33$ with $s_d = .168$ demarcating $S_*$ and $S^*$. In both panels of Figure 6, $I^*$ and $I_*$ have the boundaries given by Theorem 4 in one simulation instance. See also Figure S1.

Notice that if $x^* = 1$, then $x_*$ must also be equal to one, and in that case $I\to 0$ and $S\to 0$ $a.s.$ If $x^* = 0$, then $1-s_d$ is a divider of the $I^*$ and $I_*$. Moreover, if $x^* < 1-s_d$, then $\limsup\nolimits_{t\to \infty} I(t) > 0$ and so $\liminf\nolimits_{t\to \infty} I(t) > 0$ or otherwise $\lim_{t\to \infty} I(t) = 0$ since $0$ is an absorbing state for $I$. Note that $s_d$ is always greater than or equal to $\frac{1}{R_0}$ and are equal when $\sigma_1^2 = 0$. We can infer that $1-s_d$ is a stochastic herd immunity threshold $HIT^s$. Interestingly, $HIT^s$ is less than or equal to the deterministic herd immunity threshold $1-\frac{1}{R_0}$.

Note that simulating an almost sure result requires a large number of simulation runs in which the proportion of those simulation runs revealing the result is almost one. Also, we estimated the limit supremum and infimum by long-term simulations in which the cumulative maximum and minimum values of the time-reverse of the process are found to become stable. We notice that they reach stability over a period of time but lose that stability in the last period of the simulation run, since the maximum and minimum are found over a small period of time at the end of the simulation.

The following theorem shows that the limiting temporal variability for fluctuations around levels less than or equal to the equilibria is affected by the noise in the stochastic replicator dynamics representing the parental behavior.

Theorem 5. Let $(S_e, I_e, x_e)$ be an equilibrium of $\mathcal{E}_5$ type in which $x_e\neq 0, 1$. If $\frac{1}{2}\frac{2 \mu + \gamma}{\beta} \sigma_1^2 I_e < \mu$ and $\delta < \frac{1}{4} k \sigma^2$, where $\sigma^{2} = \sigma_{2}^{2}+\sigma_{3}^{2}$, then

almost surely, where $m = \min(\mu -\frac{1}{2}\frac{2 \mu + \gamma}{\beta} \sigma_1^2 I_e, \, \mu +\gamma, \, \mu \left(\frac{1}{2} k \sigma^2 - 2 \delta \right))$.

Proof. Let $a = \kappa \sigma_{3}^{2} -\delta - \omega$ and $b = \kappa \sigma^2 -2 \delta$. Notice that we can write $a = bx_e-I_e$ and $I + a - bx = (I - I_e) - b(x-x_e)$ since $I_e, x_e$ are zeros of the drift in the stochastic replicator Eq. Consider the following three non-negative functions:

(Note: $V_1$ is the relative entropy.) Thus, Itô formula (see Lemma S2 in the Supplementary) implies

Moreover,

but since $\mu = \mu x_e + \mu S_e + (\mu + \gamma) I_e$ where $S_e, I_e, x_e$ are zeros of drift of the first Eq in (2.5), then

Finally,

Now, let us introduce the non-negative function

Thus,

But, since $- \mu (S - S_e)(x -x_e)\leq \mu S_e x + \mu S x_e\leq \mu S_e + \mu x_e$, while

then

Now, if $dV\leq (f(S_e, I_e, x_e)-g(S, I, x))dt+\sigma(S, I, x)dW_t$ for some positive functions $V, f, g$, then

almost surely. But then

Thus,

almost surely and the result follows.

If $\sigma_1^2 = 0$, as in Figure 6(a), then

almost surely.

4.   Discussion

In this paper, we introduced a novel model of stochastic replicator dynamics (RD), from stochastic game theory, to depict human choice with noise in their perceived utility function. Including white noise in the utility function serves as a way to model irrationality and stochastic deviations from the classical rational utility. Introducing the white noise to the perceived utilities of both groups gives rise to a mutation term in the drift (deterministic) component in the stochastic RD. The mutation term depicts irrational human choices in which they explore other, possibly non-optimal, behaviors. In addition, the stochastic component of the RD reflects the unpredictable nature of human behavior and the various biases and constraints that influence decision-making. Thus, adding stochastic deviations from the classical rational utility aligns more closely with the bounded rationality framework. To our knowledge, this is the first time that parental choice of vaccination for their children has been modeled using stochastic replicator dynamics to describe their bounded rationality. See also ^[27,44,45] for other deterministic models of bounded rationality.

We also included noise in disease transmission to investigate its effect on vaccination decision-making. Although the latter has been investigated in the literature, this is the first time that the coupled system of stochastic differential Eqs of SIR and RD has been studied. The coupling takes place as the stochastic process of the proportion of infected children appears in the utility functions and as the stochastic process of the vaccine uptake appears in the first Eq of the stochastic SIR. The model shows new dynamics that generalize to those investigated in ^[28].

The presence of stochastic noise in the perceived utilities has no effect on the disease transmission if the stochastic basic reproduction number is less than one and the noise in the disease transmission is small enough. But when those conditions change, stochastic noise in perceived utility can affect the dynamical behavior of the disease due to the difference in the magnitude of that noise in the vaccinator and non-vaccinator parents in a way that can counter peer pressure. The larger the magnitude of that noise in the vaccintor's perceived utility, the more the likelihood that parents will not vaccinate their children, and vice versa for non-vaccinators. A clear and consistent view of the utility of vaccinators is required to achieve disease eradication when the disease has a large stochastic basic reproduction number. Disease eradication becomes uncertain when both magnitudes of noises approach the same value. It also depends on the initial acceptance of the vaccine.

We noticed here the emergence of new stochastic herd immunity thresholds that are less than the deterministic one. If vaccination uptake is lower than those new thresholds, then the disease will persist.

The effect of stochasticity on vaccine acceptance and disease eradication could be controlled through a stochastic optimal control. A discount in the cost of vaccination, including but not limited to ease of access to a safe vaccination as well as lowering its price, can achieve full coverage of the vaccine. That is true even when epidemiological, social, and cognitive parameters induce disease persistence. It is shown here that an optimal discount is not required for the full period and could be gradually eased up. The optimal control policy takes into account that vaccination efforts have practical and resource-related bounds such as finite resources and supply.

Although the underlying model attempts to capture the population-level effects of shared or systematic uncertainties, it does not account for heterogeneity in the individual's perceptions of payoffs. The assumption of shared or systematic stochasticity aligns with situations in which systemic external influences, e.g., environmental, cultural, informational, and economical, dominate the individual differences, especially in tightly connected or homogeneous populations. It may be less applicable to diverse populations with significant individual variability, which emphasizes the need for models that balance collective and individual-level influences on disease spread. Future work in stochastic modeling of decision-making in behavioral epidemiology should extend itself to incorporating individual-specific (occasional) stochasticity, which would allow for the exploration of diverse behaviors and richer dynamics at both individual and population levels.

In our simulations, where time is measured in years, the proportions of the susceptible, infected, and vaccine acceptors sometimes take a long time to converge to equilibrium under certain parameter settings. This slow convergence reflects the slow mixing dynamics of the system. In future work, it is warranted to find or estimate parameter values that reflect the empirically observed behavior of these variables in real life.

5.   Conclusions

Adopting stochastic game theory in epidemiological modeling assists in understanding how diseases spread and in forecasting and managing public responses to health crises. This approach provides a more nuanced view of vaccine uptake through human behavior and its impact on public health, establishing behavior management as an essential element of contemporary disease control strategies.

Acknowledgments

The authors express their sincere gratitude to the editor and the two referees for their invaluable insights and constructive suggestions, which significantly improved the quality and clarity of this manuscript.

Use of AI tools declaration

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

Conflict of interest

The authors declare no conflict of interest.

S1.   Supplementary

In the Supplementary, we will present some important definitions, lemmas, and propositions, some of which can be found in the literature. The rest are proven here.

S1.1.   Definitions and supporting lemmas

Lemma S1. The generalized L'Hôpital rule ^[46]

Let $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq 0}, \mathbb{P})$ be a complete probability space with filtration $\{\mathcal{F}_t\}_{t\geq 0}$ with respect to which all of the stochastic processes introduced in the following are defined.

Lemma S2 (Itô's formula). Let the one-dimensional stochastic differential Eq

If $f\left(t, x\right)$ is a twice-differentiable scalar function, then

A d-dimensional stochastic differential Eq

where the function $F(t, x) = \left(f_{i}(t, x)\right)_{1\leq i\leq d}$ is a $d\times 1$ vector defined on $[0, \infty]\times \mathbb{R}^d$ and $G(t, x) = \left(g_{i, j}(t, x)\right)_{1\leq i\leq d, 1\leq j\leq n}$ is a $d\times n$ vector defined on $[0, \infty]\times \mathbb{R}^d$ are both locally Lipschitz functions in $x$. Let $W(t)$ be an $n-$ dimensional Brownian motion. Let $V$ be a function in $C^{1, 2}([0, \infty]\times \mathbb{R}^d, [0, \infty))$ and let the differential operator $L$ be defined by

where $DV(t, x) = \dfrac{\partial V(t, x)}{\partial x}$ is the gradient of $V(t, x)$ and $H(t, x) = \dfrac{\partial^2 V(t, x)}{\partial x^2}$ is the Hessian matrix of $V(t, x)$. Itô's formula for $X(t)$ and such function $V$ states that

for $t\geq 0$.

The following lemma can be found in ^[40].

Lemma S3. Let $\{M_t:t\geq 0\}$ be a real-valued continuous local martingale with respect to $\{\mathcal{F}_t\}_{t\geq 0}$ such that $M_0 = 0$. If $\limsup\nolimits_{t \to \infty} \frac{\mathbb{E}(M_t^2)}{t} < \infty$ then $\lim_{t \to \infty} \frac{M_t}{t} = 0$ $a.s.$

The following lemma is an extension to the lemmas that were introduced in the Supplementary of ^[47]. It is new to our knowledge.

Lemma S4. Consider a function $f\in C([0, \infty)\times \Omega, (0, \infty))$. Let $a(t)$ be a sequence such that $\lim_{t \to \infty} \dfrac{a(t)}{t} = a$ a.s. and let $a$ and $b$ be two positive constants such that

for all $t\geq 0$ a.s. Then

Proof. We will prove the "$\geq$" part and the "$\leq$" follows using the same approach. Let $g(t): = b\, f(t)$ and so the inequality becomes

Let $\lim_{t \to \infty} \dfrac{a(t)}{t} = a$ a.s. hence for arbitrary $\epsilon > 0$ but small enough so that there exists $\Omega_\epsilon\in \mathcal{F}$ with $\mathbb{P}(\Omega_\epsilon)\geq 1-\epsilon$ and $0 < a-\epsilon \leq \dfrac{a(t)}{t}(\omega)\leq a+\epsilon$ for all $t\geq t_0(\omega)$ for each $\omega \in \Omega_\epsilon$. Let $G(t) = \int_0^t g(s) ds$ for all $t\geq 0$, then $dG(t) = g(t)dt$ and

for all $t\geq 0$. Therefore,

and $e^{a(t)}\geq e^{(a-\epsilon)t}$ on $\Omega_\epsilon$. Thus, for $t > t_0$

Hence, by L'Hôpital rule (see Lemma S1)

and so

for arbitrary $\epsilon > 0$, thus

S1.2.   Existence, uniqueness, and boundedness

Here, we will prove the positive invariance of the solution set, $\mathcal{S} = \{(s, i, x)\in \mathbb{R}_+^3: 0\leq s+i \leq 1 \text{ and } 0\leq x \leq 1\}$. It is not straightforward since the drift and volatility in the stochastic Eqs are locally Lipschitz continuous but not growing linearly. That might cause a blow-up in the solutions of the stochastic differential Eqs.

Lemma S5. If model (2.5) has initial state $(S(0), I(0), x(0))\in \mathcal{S}$, then the model has a unique solution $(S(t), I(t), x(t))\in \mathcal{S}$ for all $t\geq 0$ with probability one.

Proof. All of the relationships in this theorem are true point-wise. By construction of the stochastic replicator Eq $0\leq x(t)\leq 1$ for all $t > 0$ if $0\leq x(0)\leq 1$. Also, $d(S+I+R) = 0$, thus if $0\leq S(0)+I(0)\leq 1$ and

for all $t > 0$. Thus, for any $t > 0$,

On the other hand,

for all $t > 0$. Thus, for any $t > 0$,

The processes $S$, $I$, and $x$ are positive bounded processes. Moreover, they have bounded limits are shown in the following proposition.

Proposition S1. If the initial state of the model (2.5) is such that $(S(0), I(0), x(0))\in \mathcal{S}$, then the following inequality always holds true.

Proof. Since

then

Inequality (S2) follows from using Lemma S5 and the generalized L'Hôpital rule (see inequality (S1) in the Supplementary) applied to the terms $e^{-\mu t}\int_{0}^{t} e^{\mu s}x(s)ds$ and $e^{-\mu t}\int_{0}^{t} e^{\mu s}I(s)ds$.

S1.3.   Stochastic optimal control problem to achieve vaccine acceptance

For convenience, we define the vector $y(t) = [S(t), I(t), x(t)]^\top$. The system can be re-written as the SDE vector form

with initial condition $y(0) = [S(0), I(0), x(0)]^\top = y_0$ where the functions $f$ and $g$ are vectors with the following components

and

respectively. We consider the quadratic cost functional:

for some constants $\alpha_1, \alpha_2, \alpha_3 \geq 0$. The goal is to find an optimal control $u^*(t)$ such that $J(u) \leq J(u^*)$, for all $u \in U$, where the set of admissible controls is $U = \{u(t): u(t) \in [0, u_{\max}], \text{ for all } t \in (0, T_f] \}$ for some $0 \leq u_{\max} < 1$, and $T_f$ is the termination time.

To use the stochastic maximum principle, we define Hamiltonian $H(y, u, p, q)$ by

where $\langle \mathbf{u}, \mathbf{v} \rangle = \sum_{\mathbf{i} = \mathbf{1}}^\mathbf{n} \mathbf{u}_\mathbf{i} \mathbf{v}_\mathbf{i}$ is the Euclidean inner product and $l(y, u) = \alpha_1 S + \alpha_2 I + \frac{\alpha_3}{2} u^2$, while $p = [p_1, p_2, p_3]^\top$ is the adjoint vector of costates and $q = [q_1, q_2, q_3]^\top$ are costates volatilities. Hence, the associated Hamiltonian is given by

According to the stochastic maximum principle,

where the optimal states satisfy

with the backward stochastic differential Eqs of the costates

where $u^*$ solves

and exists since $H$ is concave in $u$. Thus, the optimal control is given by

The function $p_3(t)$ is found by solving the system of backward stochastic differential Eqs (BSDEs) given by

such that $p_i(T_f) = 0$ for $i = 1, 2, 3$. While $x^*(t)$ solves Eq (3.10) at $u^*(t)$ given by (S6). Including the reduction of susceptible in the objective function results in a longer reduction in the vaccination cost than when excluded from the objective.

S1.4.   More simulations