In this paper, we model dynamics of pediatric vaccination as an imitation game, in which the rate of switching of vaccination strategies is proportional to perceived payoff gain that consists of the difference between perceived risk of infection and perceived risk of vaccine side effects. To account for the fact that vaccine side effects may affect people’s perceptions of vaccine safety for some period of time, we use a delay distribution to represent how memory of past side effects influences current perception of risk. We find disease-free, pure vaccinator, and endemic equilibria and obtain conditions for their stability in terms of system parameters and characteristics of a delay distribution. Numerical bifurcation analysis illustrates how stability of the endemic steady state varies with the imitation rate and the mean time delay, and this shows that it is not just the mean duration of memory of past side effects, but also the actual distribution that determines whether disease will be maintained in the population at some steady level, or if sustained periodic oscillations around this steady state will be observed. Numerical simulations illustrate a comparison of the dynamics for different mean delays and different distributions, and they show that even when periodic solutions are observed, there are differences in their amplitude and period for different distributions. We also investigate the effect of constant public health information campaigns on vaccination dynamics. The analysis suggests that the introduction of such campaigns acts as a stabilizing factor for endemic equilibrium, allowing it to remain stable for larger values of mean time delays.

Control of infectious diseases using vaccination and other public health methods presents a fascinating yet challenging problem from the perspective of understanding and modeling collective behavior. People are known to modify their attitudes to vaccines due to personal experience, information spread through word-of-mouth, as well as messages from various public health bodies, and this can have a profound impact on their vaccination choices. There have been several notable vaccine scares in the last 30 years, with many countries still suffering from measles and mumps outbreaks due to insufficient coverage with measles, mumps, and rubella (MMR) vaccines. More recently, misinformation about COVID-19 vaccines led to significant challenges in their adoption in many countries despite excellent efficacy and safety profiles. This paper will use a known framework of evolutionary game theory and, more specifically, imitation games, where the change in the proportion of vaccinators increases with their perceived payoff that consists of the difference between their risk of infection and risk of vaccine side effects. We pay particular attention to the fact that vaccine side effects may appear not immediately after vaccination, and their perceived risk may also not dissipate instantly, which is modeled using distributed time delays. We identify disease-free, endemic, and other steady states of the model and study their stability depending on model parameters and delay distributions. Whereas stability of the disease-free steady state is distribution-independent, for the endemic equilibrium, we find that it is not only the mean time delay, but also the actual shape of the delay distribution that determines whether there will be some steady level of infection or periodic oscillations in disease incidence. We also explore the effect of public health campaigns on controlling the disease, as well as on stability of endemic equilibrium. Numerical simulations are performed to illustrate the behavior of the model in different dynamical regimes.

In the context of understanding the dynamics and control of infectious diseases, a significant role is played by human behavior that can have major effects on epidemic dynamics. This includes people reducing their contacts and otherwise modifying their behavior, once they become aware of the ongoing outbreak, as well as using vaccines to prevent infection. At the same time, one should note that information and awareness can result not only in the positive outcome in terms of reducing spread of disease, but it can also have a major detrimental effect when it causes panic and/or undermines public confidence. One example of the latter is the uncontrolled spread of plague during the 1994 outbreak in one of the states in India, where people fled the endemic and, as a result, carried the disease with them and infected other parts of the country.1 Another notable example of vaccine fears with global ramifications is the case of a 1998 paper in Lancet alleging a possible connection between MMR vaccines and autism. Despite this paper being subsequently retracted, and the claim itself being thoroughly debunked with tens, if not hundreds, of controlled studies, and hundreds of thousands of people, even now, there are significant numbers of young people who are not vaccinated against measles, mumps, and rubella in different countries. Another more recent example concerns COVID-19 pandemic and vaccination against SARS-CoV-2, where misinformation and some rather wild conspiracy theories led to a notable negative impact on vaccine take-up.2–4 Besides fears of vaccine-induced side effects, an important role is also played by the public perception and the way immunization campaign is portrayed in public media. For example, a negative press coverage of an HPV vaccination campaign in Romania, which claimed that vaccination would encourage promiscuous behavior among young girls,5 resulted in a very low take-up of vaccination, leading to Romania having the highest rates of cervical cancer in Europe.6 

When discussing vaccination, important concepts are those of herd immunity and free-riding. Herd immunity refers to the idea that if a sufficiently high percentage of population is vaccinated, even those who are not vaccinated are benefiting from protection, because an outbreak will not occur, even if infection is introduced in the population due to the pool of susceptible people being very small.7–9 Mathematically, if the infectious disease is characterized by the basic reproduction number R 0, then, subject to a number of simplifying assumptions, such as heterogeneity of contacts, absence of age structure, etc., the critical percentage of population that needs to be vaccinated in order to achieve herd immunity can be estimated as p c = 1 1 / R 0. Measles, which is known to be the most contagious infectious disease to date, has the value of R 0 estimated between 12 and 1810 (with some more recent estimates suggesting a range of between 6.51 and 45.4 in European countries11), with the herd immunity threshold of 94.4% for R 0 = 18. In 2015–2017, over 75% of European countries had a lower percentage of people with two doses of an MMR vaccine,12 which explains continuing outbreaks of measles in these and other countries. Free-riding refers to the strategy taken up by some people where they are relying on herd immunity for their protection, instead of actually getting vaccinated themselves, thus avoiding potential vaccine side effects (VSEs).13–16 This approach is different from simply refusing vaccination on some religious, philosophical, or other grounds.18 Several studies have shown that free-riding makes it impossible to completely eliminate infection from population.13,14,16–19

In terms of mathematical modeling of epidemic dynamics with account for changes in human behavior, a variety of different approaches have been used that focused on behavioral changes resulting from awareness spread through direct social contacts or via public health information campaigns. A monograph by Manfredi and d’Onofrio20 provides an excellent overview of many techniques used in such modeling (see also a more recent review paper by Wang et al.21). One popular approach to studying vaccination dynamics, which has been successfully used in various contexts over the last almost 20 years, is based on evolutionary game theory. This approach originated in the work of Bauch14,22 that proposed modeling a proportion of vaccinators (people willing to vaccinate) using an imitation game, where non-vaccinators imitate vaccinators and take up vaccination at a rate that is proportional to their perceived payoff from adopting a vaccinator strategy. This payoff is then taken to be the difference in the perceived risk of disease, which is usually taken to be proportional to disease incidence, or some other growing function of it, and the perceived risk of VSE, which is represented as a linear (or some other) function of vaccination coverage. In the simplest case, where risk of disease is proportional to disease incidence and the risk of side effects is proportional to vaccination coverage, dynamics of the proportion of vaccinated people is described by the well-known replicator equation from evolutionary game theory.23 d’Onofrio et al.24 included in the model additional perceived payoff from vaccination that arises due to a public health campaign, which, for simplicity, they considered to be constant. Buonomo et al.25 studied the effects of seasonality in the transmission rate in such a model, while Buonomo et al. have used optimal control theory to obtain optimal time profiles of public health awareness programs in the case of constant transmission26 or seasonally forced transmission.27 

Vaccine-induced side effects often appear not immediately, but rather some time after vaccination, and once information about the side effects becomes known, this has the potential to affect vaccine take-up for some considerable time afterward. One way to account for these effects in mathematical models is to consider the current proportion of people deciding to take up vaccination being a positive decreasing function of a variable representing the perceived risk of vaccine, which can itself be delayed.19,28 Another approach, which was briefly mentioned in d’Onofrio et al.,17 uses the imitation game model14 but allows current perceived risk of vaccination to depend on past vaccine side effects, included with a delay kernel. While the main bulk of that work focused on an unlagged model, numerical simulations with exponentially fading memory demonstrated the possibility of sustained oscillations around the endemic equilibrium. Following similar methodology, in this paper, we consider a pediatric infectious disease, such as measles or mumps, which is controlled by a vaccine that is assumed to provide 100% effectiveness and a lifelong immunity. The disease itself is also assumed to confer lifelong immunity, which means that disease dynamics can be modeled using an SIR-type model. Vaccine is administered to newborns only, and the parents make a choice of whether or not to vaccinate depending on available information about vaccine risks vs risks of infection, and in doing so, parents imitate other parents who have vaccinated their children. The dynamics of such a model is then given by the following system of equations:
S ˙ = μ ( 1 p ) β S I μ S , I ˙ = β S I ( μ + v ) I , p ˙ = k p ( 1 p ) Δ E ,
(1)
where S and I are proportions of susceptible and infected individuals; μ = 1 / L is the population birth rate taken to be the same as the death rate, where L is the life expectancy at birth; β is the disease transmission rate; v is the recovery rate; and k is the imitation rate, at which individuals sample and adopt an alternative vaccination strategy depending on the payoff gain. Δ E = P v P i is the difference in perceived payoffs for vaccinated and unvaccinated individuals, where P v = r v and P i = r i, with r v and r i being the perceived risk of VSE and perceived risk of infection, respectively. Similarly to other models,14,17,24 we will assume perceived risk of infection r i to be proportional to disease prevalence; i.e., r i = ω I. Perceived risk of vaccination r v is taken in the form that depends on past vaccine-related side effects,
r v ( t ) = α 0 g ( s ) p ( t s ) d s = α t g ( t z ) p ( z ) d z ,
where g ( s ) is the delay kernel, which quantifies how memory of past VSE is influencing current perception of risk. The delay kernel g ( s ) is assumed to be positive-definite, i.e., g ( s ) 0 for s 0, and normalized to unity, i.e., 0 g ( s ) d s = 1. The complete model now has the form
S ˙ = μ ( 1 p ) β S I μ S , I ˙ = β S I ( μ + v ) I , p ˙ = k p ( 1 p ) [ I α 0 g ( s ) p ( t s ) d s ] ,
(2)
where we have introduced rescaled parameters
k ¯ = k ω , α ¯ = α ω
and removed bars for simplicity.
Specific delay kernels we will consider below include a discrete delay distribution g ( s ) = δ ( s τ 0 ) that assumes that the perceived risk of VSE at time t is determined by VSE that occurred some fixed period of time τ 0 ago, i.e., at time t τ 0; Gamma distribution g ( s ) = s n 1 σ n e σ s / ( n 1 )!, which for n = 1 (weak kernel, exponentially fading memory17) and n = 2 (strong kernel) turns into g ( s ) = σ e σ s and g ( s ) = σ 2 s e σ s, respectively, and an acquisition-fading (AF) kernel28  g ( s ) = ( e s / T 1 e s / T 2 ) / ( T 1 T 2 ). In the limit T 1 0, T 2 = T (or T 2 0, T 1 = T), this kernel transforms into the weak kernel with σ = 1 / T, and in the limit T 1 T 2 = T, AF kernel becomes the strong kernel with the same σ = 1 / T. These delay kernels are illustrated in Fig. 1, with the exception of discrete delay kernel given by the Dirac δ-function at τ 0. The mean time delay for each of these delay distributions can be found as
τ = 0 s g ( s ) d s ,
and it gives τ = τ 0 for the case of discrete delay, τ = n / σ for a Gamma-distributed delay kernel, and τ = T 1 + T 2 for the AF kernel.
FIG. 1.

Weak, strong, and acquisition-fading distribution kernels with the mean time delay τ = 2.

FIG. 1.

Weak, strong, and acquisition-fading distribution kernels with the mean time delay τ = 2.

Close modal

The work by d’Onofrio et al.17 mentions the model (2) with the exponentially fading memory (weak Gamma distribution), and while no analytical results were obtained, for this specific delayed kernel, the authors presented a bifurcation diagram for the endemic steady state in the plane of inverse mean time delay and the imitation rate k. They also illustrated regimes where endemic equilibrium is stable or unstable by showing the corresponding numerical solutions. Unfortunately, the highlighted parts of their bifurcation diagram (Fig. 6 in Ref. 17) are referred to in the caption to that figure as stability regions, and in the accompanying text as instability regions, and without analytical results supporting this, the reader is left confused. To clarify this confusion, in this paper, we derive some general kernel-independent results and perform both analytical and numerical investigations to analyze stability of all steady states of the model (2) not only for the exponentially fading memory, but for all four of the above-mentioned delay kernels.

The remainder of this article is organized as follows. In Sec. II, we identify steady states of model (2) and derive conditions for their biological feasibility and stability. Section III contains the results of numerical bifurcation analysis, where we explore numerically stability of endemic equilibrium, and we solve the model numerically, which allows us to compare and contrast behavior of the system for the same values of parameters but different delay distributions. In Sec. IV, we extend the model to include the impact of public health information campaign and analyze its impact on disease control and on stability properties of an endemic steady state. The article concludes in Sec. V with a discussion of results.

For any values of parameters, model (2) has steady states E 0 = ( 1 , 0 , 0 ) and E 1 = ( 0 , 0 , 1 ). It can also have steady states
E 2 = ( μ + v β , μ [ β ( μ + v ) ] β ( μ + v ) , 0 )
and
E 3 = ( S , I , p ) = ( μ + v β , μ α [ β ( μ + v ) ] β [ μ + α ( μ + v ) ] , μ [ β ( μ + v ) ] β [ μ + α ( μ + v ) ] ) .
The steady state E 0 is the disease-free steady state of the model, E 1 is the so-called pure vaccinator equilibrium,14,24,27 E 2 is the endemic state in the absence of vaccination, and E 3 is the fully endemic steady state with vaccination. Both steady states E 2 and E 3 are only biologically feasible if β > ( μ + v ).
Stability of any steady state E ^ ( S ^ , I ^ , p ^ ) of the model (2) is determined by the roots λ of the characteristic equation
| μ β I ^ λ β S ^ μ β I ^ β S ^ ( μ + v ) λ 0 0 k p ^ ( 1 p ^ ) k 1 k 2 { L g } ( λ ) λ | = 0 ,
(3)
where
k 1 = k ( 1 2 p ^ ) ( I ^ α p ^ ) , k 2 = k α p ^ ( 1 p ^ ) ,
and
{ L g } ( λ ) = 0 g ( s ) e λ s d s
is the Laplace transform of the kernel g ( s ). For discrete delay, we have
{ L g } ( λ ) = e λ τ 0 ,
and, similarly, for the Gamma distribution,
{ L g } ( λ ) = ( σ λ + σ ) p ,
and for the AF kernel,
{ L g } ( λ ) = 1 ( 1 + λ T 1 ) ( 1 + λ T 2 ) .
For the disease-free steady state E 0, it immediately follows from (3) that this steady state is stable, provided
R 0 = β μ + v < 1 ,
(4)
unstable if R 0 > 1, and undergoes a steady-state bifurcation at R 0 = 1. Here, R 0 is the basic reproduction number, which, importantly, only depends on disease parameters but is independent of delay distributions or parameters k and α characterizing vaccination. Boundary equilibrium E 1 has eigenvalues λ 1 = μ, λ 2 = ( μ + v ), and λ 3 = α k, the last of which is positive, implying that this steady state is always unstable. Similarly, one of the eigenvalues of the steady state E 2 is given by λ = k I , which indicates that as long as this steady state exists, it will also always be unstable. We note that the condition of instability of disease-free steady state R 0 > 1 is also the condition that ensures biological feasibility of the steady states E 2 and E 3.
The characteristic equation for the endemic equilibrium E 3 has the form
k μ β p ( 1 p ) I + [ k α p ( 1 p ) { L g } ( λ ) + λ ] × [ λ 2 + λ ( μ + β I ) + β 2 S I ] = 0 ,
(5)
and for very small imitation rates 0 < k 1, this reduces to a simple quadratic equation λ 2 + λ ( μ + β I ) + β 2 S I = 0, whose both roots have negative real parts. This suggests that for a sufficiently small k, the steady state E 3 is stable for all distribution kernels independently of the mean time delay.
For discrete delay, characteristic Eq. (5) transforms into a transcendental equation
P ( λ ) + Q ( λ ) e λ τ 0 = 0 ,
(6)
where P ( λ ) and Q ( λ ) are, respectively, cubic and quadratic polynomials of the form
P ( λ ) = λ 3 + λ 2 ( μ + β I ) + β 2 S I λ + k μ β p ( 1 p ) I
and
Q ( λ ) = k α p ( 1 p ) [ λ 2 + λ ( μ + β I ) + β 2 S I ] .
At τ 0 = 0, a transcendental characteristic equation (6) reduces to a cubic
λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0 ,
with coefficients
a 2 = μ + β I + k α p ( 1 p ) , a 1 = β 2 S I + k α p ( 1 p ) ( μ + β I ) , a 0 = k β p ( 1 p ) I ( μ + α β S ) .
Since all coefficients a 0, a 1, and a 2 are positive, the condition for stability of the endemic steady state E at τ 0 = 0 can be found from the Routh–Hurwitz criterion as
k α p ( 1 p ) [ ( β I + μ ) [ k α ( 1 p ) + 1 ] β μ p ] + β 2 S I ( β I + μ ) > 0.
(7)
It is easy to check that λ = 0 is not a solution of the characteristic equation (6), and hence, if the steady state E is stable for τ 0 = 0, for τ 0 > 0, its stability can only change through a Hopf bifurcation, where a pair of complex conjugate eigenvalues crosses the imaginary axis with a non-zero speed. Looking for solutions of the characteristic equation (6) in the form λ = i ζ, substituting this into the equation, and separating real and imaginary parts, we obtain an equation for the Hopf frequency ζ,
η 3 + [ ( β I + μ ) 2 2 β 2 S I ( k α p ( 1 p ) ) 2 ] η 2 + [ ( β 2 S I ) 2 2 k μ β p ( 1 p ) I ( β I + μ ) ( k α p ( 1 p ) ) 2 ( ( β I + μ ) 2 2 β 2 S I ) ] η + [ k β p ( 1 p ) I ] 2 [ μ 2 α 2 β 2 ( S ) 2 ] = 0 ,
where η = ζ 2. If the free term in this equation is negative, which happens whenever the following condition holds:
α β S > μ α v + μ ( α 1 ) > 0 ,
there will exist at least one real positive root η 0 of this equation, which would give the corresponding Hopf frequency as ζ 0 = η 0. The critical value of the time delay τ 0, at which this Hopf bifurcation would occur, can then be found as
τ 0 c = arctan k ^ ( ζ 0 2 β 2 S I ) 2 + ζ 0 ( μ + β I ) z 0 ζ 0 2 ( μ + β I ) k ^ ( ζ 0 2 β 2 S I ) z 0 ,
where
k ^ = k α ζ 0 p ( 1 p ) , z 0 = k μ β p ( 1 p ) I ζ 0 2 ( μ + β I ) .
For a Gamma-distributed kernel, the characteristic equation (5) transforms into a polynomial equation of degree ( n + 3 ),
k μ β p ( 1 p ) I ( λ + σ ) n + [ k α p ( 1 p ) σ n + λ ( λ + σ ) n ] × [ λ 2 + λ ( μ + β I ) + β 2 S I ] = 0.
(8)
Another way of arriving at this equation is to use a linear chain trick to transform the original model (2) with a distributed delay into a system of coupled ordinary differential equations. To illustrate this process for the weak kernel n = 1, we introduce an auxiliary variable M ( t ) as
M ( t ) = 0 σ e σ s p ( t s ) d s .
Then, the system (2) with a weak kernel can be equivalently rewritten as
S ˙ = μ ( 1 p ) β S I μ S , I ˙ = β S I ( μ + v ) I , p ˙ = k p ( 1 p ) [ I α M ] , M ˙ = σ p σ M .
(9)
The endemic steady state E 3 transforms into the steady state E ~ 3 = ( S , I , p , p ) of the extended system (9). Jacobian of linearization around this steady state has the form
J E ~ 3 = ( μ β I β S μ 0 β I 0 0 0 0 k p ( 1 p ) 0 k α p ( 1 p ) 0 0 σ σ ) ,
and the characteristic equation det [ J E ~ 3 λ I ] = 0 has exactly the same form as (8). Because this equation is polynomial, we can again resort to Routh–Hurwitz criteria to obtain explicit conditions for stability of endemic equilibrium E 3. Rewriting Eq. (8) with n = 1 in the form
λ 4 + b 3 λ 3 + b 2 λ 2 + b 1 λ + b 0 = 0 ,
with coefficients
b 3 = ( β I + μ + σ ) , b 0 = k β σ p ( 1 p ) I ( α β S + μ ) , b 2 = k α σ p ( 1 p ) + σ ( β I + μ ) + β 2 S I , b 1 = k p ( 1 p ) [ α σ ( β I + μ ) + μ β I ] + σ β 2 S I ,
and since these coefficients are always positive, conditions for stability of E 3 are given by
b 3 b 2 b 1 > 0 , ( b 3 b 2 b 1 ) b 1 b 3 2 b 0 > 0.
(10)
Applying a linear chain trick to the cases of strong and AF kernels results in the following systems of ODEs:
{ S ˙ = μ ( 1 p ) β S I μ S , I ˙ = β S I ( μ + v ) I , p ˙ = k p ( 1 p ) [ I α M 2 ] , M 1 ˙ = σ p σ M 1 , M 2 ˙ = σ M 1 σ M 2 , { S ˙ = μ ( 1 p ) β S I μ S , I ˙ = β S I ( μ + v ) I , p ˙ = k p ( 1 p ) [ I α M 2 ] , M 1 ˙ = a 1 p a 1 M 1 , M 2 ˙ = a 2 M 1 a 2 M 2 ,
(11)
where a 1 , 2 = 1 / T 1 , 2, and we again note that setting a 1 = a 2 = σ transforms the AF kernel into a strong kernel. The characteristic equation for the AF kernel is a quintic
[ k α p ( 1 p ) + λ ( 1 + λ T 1 ) ( 1 + λ T 2 ) ] × [ λ 2 + λ ( μ + β I ) + β 2 S I ] + k μ β p ( 1 p ) I ( 1 + λ T 1 ) ( 1 + λ T 2 ) = 0 ,
and its stability can also be investigated using Routh–Hurwitz criteria, though expressions would look more cumbersome.

To explore how stability of the endemic steady state E 3 changes with parameters, we fix the values of μ = 1 / ( 70 365 ) days 1, which corresponds to an average life expectancy of 70 years, and v = 1 / 7 days 1, which is representative of measles infection. We also fix the value of basic reproduction number R 0 = 10, which then determines the disease transmission rate β through β = R 0 ( μ + v ). Figure 2 illustrates how stability of endemic steady state E 3 varies with the mean time delay τ and the imitation rate k for the four delay distributions we consider. For all delay distributions, we note a small region characterized by very low values of k, where the endemic steady state is stable for any values of mean delay, in agreement with our discussion of the characteristic equation (5).

FIG. 2.

Regions of stability and instability of the endemic steady state E 3 for different delay distribution kernels depending on the mean time delay τ and imitation rate k. (a) Discrete delay kernel. (b) Weak Gamma distribution. (c) Strong Gamma distribution. (d) AF kernel with T 1 , 2 = ( 0.5 ± 0.4 ) τ. The color code denotes max[Re( λ)], with stable regions shown in blue-green and unstable shown in yellow-red. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, α = 0.002.

FIG. 2.

Regions of stability and instability of the endemic steady state E 3 for different delay distribution kernels depending on the mean time delay τ and imitation rate k. (a) Discrete delay kernel. (b) Weak Gamma distribution. (c) Strong Gamma distribution. (d) AF kernel with T 1 , 2 = ( 0.5 ± 0.4 ) τ. The color code denotes max[Re( λ)], with stable regions shown in blue-green and unstable shown in yellow-red. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, α = 0.002.

Close modal

The first thing to note is that the stability region in the case of a weak Gamma distribution, which is often assumed in epidemic models and represents exponentially decaying memory of VSE, has a structure that is very different from that for other delay distributions in that once an endemic steady state has stabilized for sufficiently high imitation rate τ, stability is never lost again. In contrast, for other delay distributions we consider, for not too large time delays, there is a range of k values, for which endemic equilibrium is stable, but for smaller and higher values of k, it is unstable, with instability occurring as a result of a supercritical Hopf bifurcation. Epidemiologically, this means that with the exception of exponentially fading memory, increasing the rate of imitation, i.e., encouraging unvaccinators to more stringently follow their vaccinator contacts, results in destabilizing an endemic steady state, leading to the emergence of periodic solutions around this steady state, some of which are illustrated in Fig. 4. Since the AF kernel reduces to a strong Gamma kernel when T 1 = T 2, to make the distinction between these two kernels clearer, we consider the situation, where T 1 , 2 = ( 0.5 ± 0.4 ) τ. The comparison of stability regions for discrete delay, a strong Gamma distribution, and such an AF kernel shows that the stability region in the plane of imitation rate k and the mean time delay τ is increased for a strong Gamma distribution and increased further still for the AF kernel. This shows that for the same rate of imitation, a stable endemic steady state, characterized by some steady level of infection, would be maintained by larger values of τ, in other words, for a longer period of memory of VSE. For all delay distributions, regardless of the value of the imitation rate, there is some minimum value of mean time delay τ such that for a higher value of τ, an endemic steady state is unstable. Notably, this minimum time delay required for destabilization of endemic equilibrium is significantly higher for exponentially fading than for other delay distributions.

FIG. 3.

(a) Stability of endemic equilibrium E 3 with an AF kernel. The color code denotes max[Re( λ)], with stable regions shown in blue-green and unstable shown in yellow-red. The black line corresponds to a mean time delay τ = T 1 + T 2 = 150. (b) Stability boundaries of endemic steady state E 3 depending on the mean delay τ and the half-difference between T 1 and T 2, with T 1 , 2 = ( τ ± W ) / 2. E 3 is stable to the left of the boundary and unstable to the right of it. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, k = 250. (a) α = 0.002. (b) α = 0.001 (blue), α = 0.004 (red), and α = 0.02 (green).

FIG. 3.

(a) Stability of endemic equilibrium E 3 with an AF kernel. The color code denotes max[Re( λ)], with stable regions shown in blue-green and unstable shown in yellow-red. The black line corresponds to a mean time delay τ = T 1 + T 2 = 150. (b) Stability boundaries of endemic steady state E 3 depending on the mean delay τ and the half-difference between T 1 and T 2, with T 1 , 2 = ( τ ± W ) / 2. E 3 is stable to the left of the boundary and unstable to the right of it. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, k = 250. (a) α = 0.002. (b) α = 0.001 (blue), α = 0.004 (red), and α = 0.02 (green).

Close modal
FIG. 4.

Numerical solution of the model (2) for different delay distribution kernels and mean time delays. Top row (a)–(d): discrete delay; second row (e)–(h): weak Gamma kernel; third row (i)–(l): strong Gamma kernel; and fourth row (m)–(p): AF kernel. First column: τ = 20 days, second column: τ = 50 days, third column: τ = 100 days, and fourth column: τ = 150 days. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, α = 0.002, k = 400.

FIG. 4.

Numerical solution of the model (2) for different delay distribution kernels and mean time delays. Top row (a)–(d): discrete delay; second row (e)–(h): weak Gamma kernel; third row (i)–(l): strong Gamma kernel; and fourth row (m)–(p): AF kernel. First column: τ = 20 days, second column: τ = 50 days, third column: τ = 100 days, and fourth column: τ = 150 days. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, α = 0.002, k = 400.

Close modal

With the AF kernel providing a family of kernels with the same mean delay τ = T 1 + T 2, in Fig. 3(a), we plot the stability region of endemic equilibrium E 3 depending on parameters T 1 and T 2 for some fixed value of imitation rate k. The diagonal T 1 = T 2 corresponds to the case of a strong Gamma distribution, and we observe that where the endemic steady state is initially unstable for a strong Gamma distribution, as the values of T 1 and T 2 deviate from the value of τ / 2, eventually, this leads to stabilization of endemic equilibrium. This is further explored in Fig. 3(b), where we write parameters T 1 and T 2 as T 1 , 2 = ( τ ± W ) / 2. Stability boundaries of endemic equilibrium in the τ W plane show that increasing the width W of the AF distribution results in increasing the range of mean time delays τ, for which endemic equilibrium is stable. Interestingly, increasing the rate α, which measures how past VSE impact current perceived vaccine risk, reduces the range of mean time delays, for which E 3 is stable, for smaller widths W of AF kernel, but increases this range of stable mean time delays for larger distribution widths W.

Figure 4 illustrates numerical solution of the model (2) for different delay distributions and different mean time delays. We observe the emergence of stable periodic solutions around the endemic steady state for increasing mean delays, which happens first for a discrete delay distribution, then for the strong Gamma distribution, and then for the AF kernel. In the case of a weak Gamma kernel, i.e., exponentially fading memory, endemic steady state remains stable for all values of the mean time delay we consider, as the Hopf bifurcation of this steady state occurs for a much higher value of τ. We note that in the cases, where the periodic solutions are observed for different delay distributions but the same values of parameters, these periodic solutions are characterized by different amplitudes and periods. This further highlights the point that it is not just the mean time delay, but rather the actual delay distribution that determines the type of dynamics that will be exhibited by the model.

In order to investigate the effect of public health awareness campaigns on take-up of vaccination and the resulting epidemic dynamics, we extend our model in a manner similar to d’Onofrio et al.24 and introduce an extra term in equation for p ˙, namely,
p ˙ = k p ( 1 p ) [ ω I α 0 g ( s ) p ( t s ) d s ] + k P Δ E p ( 1 p ) ,
where Δ E p ( t ) is the perceived payoff from public health information and k p is the rate of acceptance of public information. For simplicity, similarly to Refs. 24 and 25, we consider Δ E p to be constant and then write it in the form k P Δ E p = k γ. Rescaling parameters
k ω = k ¯ , α ¯ = α ω , γ ¯ = γ ω ,
and dropping bars, we obtain the modified model
S ˙ = μ ( 1 p ) β S I μ S , I ˙ = β S I ( μ + v ) I , p ˙ = k ( 1 p ) [ p ( I α 0 g ( s ) p ( t s ) d s ) + γ ] .
(12)
Similarly to the original model (2), this modified model can also have up to four steady states: disease-free state E ~ 0 = ( 1 γ α , 0 , γ α ), pure vaccinator state E ~ 1 = ( 0 , 0 , 1 ), and steady states
E ~ 2 , 3 = ( μ + v β , μ β x , 1 μ + v β ( x + 1 ) ) ,
with
x = 1 2 ( μ + v ) [ α ( μ + v ) + μ ] [ [ β ( μ + v ) ] [ 2 α ( μ + v ) + μ ] ± μ 2 [ β ( μ + v ) ] 2 + 4 β 2 γ ( μ + v ) [ α ( μ + v ) + μ ] ] .
We have the following connection between steady states of the modified and original models:
lim γ 0 E ~ 0 = E 0 , lim γ α E ~ 0 = E ~ 1 = E 1 , lim γ 0 E ~ 2 , 3 = E 2 , 3 .
Before looking at stability of these steady states, we note that any solution of the system (12) that starts in the region
Ω = { ( S , I , p ) R + 3 : 0 S + I 1 , 0 p 1 }
will remain in this region for all t 0, which means that Ω is a positively invariant region.
Characteristic roots of the steady state E ~ 1 are λ 1 = μ, λ 2 = ( μ + v ), and λ = k ( α γ ), suggesting that this steady state is linearly asymptotically stable for γ > α and unstable for γ < α. Next, we use the methodology developed in Buonomo et al.25 to show that for γ α, the steady state E ~ 1 is actually globally asymptotically stable. To demonstrate this, we consider the last equation of the system (12) that can be rewritten as follows:
p ˙ = k ( 1 p ) [ p ( I α 0 g ( s ) p ( t s ) d s ) + γ ] α k ( 1 p ) [ 1 0 g ( s ) p ( t s ) d s ] = α k ( 1 p ) 0 g ( s ) [ 1 p ( t s ) ] d s 0.
With p ( t ) being contained between 0 and 1, this then implies lim t p ( t ) = 1, with p monotonically increasing. Therefore, for any T > 0, there exists ε > 0 such that p ( t ) > 1 ε for t > T, and for T , we have ε 0. Considering t > T, equation for S ( t ) now takes the form
S ˙ = μ ( 1 p ) β S I μ S < μ ( ε S ) .
An equation of the form
z ˙ = μ ( ε z )
has the solution
z ( t ) = ε ( ε z 0 ) e μ t .
Since lim t z ( t ) = ε, from the comparison principle, we have S ( t ) 0 as t . In a similar way, choosing sufficiently large T, such that S ( t ) < ε = ( μ + v ) / β, we would have for t > T,
I ˙ = β S I ( μ + v ) I < [ β ε ( μ + v ) ] I < 0 ,
which then implies I ( t ) 0 as t . Taken together, this shows that for γ α, the steady state E ~ 1 is indeed globally asymptotically stable.
Linearization of the modified model (12) near the disease-free steady state E 0, which is only biologically feasible for 0 γ / α < 1, has characteristic eigenvalues λ 1 = μ, λ 2 = β ( 1 γ α ) ( μ + v ), with remaining eigenvalues satisfying the equation
λ = k [ γ α p ( 1 2 p ) α p ( 1 p ) { L g } ( λ ) ] .
(13)
Using the relation α ( p ) 2 = γ, this equation can be equivalently rewritten as
λ = k α γ ( γ α 1 ) [ 1 + { L g } ( λ ) ] = C [ 1 + { L g } ( λ ) ] ,
(14)
where
C = k α γ ( 1 γ α ) > 0.
(15)
For the case of discrete delay, the characteristic equation (14) turns into
λ = C ( 1 + e λ τ ) .
(16)
For τ = 0, the solution of this equation is λ = C < 0. Obviously, for τ > 0, there can be no real roots of this equation, as this would result in a real negative value of the left-hand side and a real positive value on the right-hand side. What remains to be checked is whether there can be complex conjugate roots of this equation. Since at τ = 0, the characteristic eigenvalue is real and negative and λ = 0 is clearly not a solution, we need to check whether a pair of complex conjugate eigenvalues can cross the imaginary axis for some τ 0 = τ 0 c. If this happens, we would have λ = ± i ζ. Substituting this into Eq. (16) and separating real and imaginary parts, we have
{ ζ = C sin ( ζ τ ) , 0 = C [ 1 + cos ( ζ τ ) ] .
From the second equation, we obtain cos ( ζ τ ) = 1, which implies ζ τ = π + 2 n π, n = 0 , ± 1 , ± 2 , , and substituting this into the first equation, we then have sin ( ζ τ ) = 0, which means that ζ = 0. Hence, we conclude that there cannot be purely imaginary roots of Eq. (16), and all of its roots have a negative real part for any values of τ.
For a Gamma-distributed delay kernel, Eq. (13) takes the form
λ = C [ 1 + ( σ λ + σ ) n ] .
In the case of a weak kernel, this equation transforms into a quadratic equation
λ 2 + ( C + σ ) λ + 2 C σ = 0 ,
whose both roots have a negative real part for any σ > 0. In the case of a strong kernel, the characteristic equation becomes a cubic
λ 3 + ( C + 2 σ ) λ 2 + σ ( σ + 2 C ) λ + 2 C σ 2 = 0.
Since all coefficients of this equation are positive, the Routh–Hurwitz criterion for stability gives the condition
σ ( 2 C 2 + 3 C σ + 2 σ 2 ) > 0 ,
which is always satisfied. Hence, as in the case of discrete delay, Eq. (13) does not have roots with a positive real part.
Finally, for the AF kernel, Eq. (13) is also cubic,
T 1 T 2 λ 3 + ( C T 1 T 2 + T 1 + T 2 ) λ 2 + [ C ( T 1 + T 2 ) + 1 ] λ + 2 C = 0 ,
and the Routh–Hurwitz criterion,
( C T 2 + 1 ) [ ( C T 2 + 1 ) T 1 + T 2 ] + T 1 > 0 ,
is always satisfied. From these calculations, we conclude that irrespective of which distribution kernel we consider, stability of the disease-free steady state E 0 is determined by the sign of the eigenvalue λ 2, which will be negative, provided
R 0 p < 1 , R 0 p = β μ + v ( 1 γ α ) = R 0 ( 1 γ α ) .
This shows that increasing the value of γ reduces the value of basic reproduction number, thus increasing stability of the disease-free equilibrium.
For endemic equilibrium, the characteristic equation has the form
k μ β p ( 1 p ) I + [ λ 2 + λ ( μ + β I ) + β 2 S I ] × [ k ( 1 p ) [ α p { L g } ( λ ) + ( I α p ) ] + λ ] = 0.
(17)
The difference from the case of no public health intervention is in the extra term ( I α p ) inside the last bracket, and this term is equal zero when γ = 0.

In Fig. 5, we have fixed the value of imitation rate k and the mean time delay τ and plotted regions of stability and instability of endemic equilibrium E ~ 3 of the modified model (12) depending on the strength γ of the public health information campaign and the rate α, at which past VSE determines current perceived vaccine risk. For the case of discrete delay, there is a very narrow region, characterized by the highest values of γ, at which an endemic steady state is still feasible, where it is stable, while in the rest of the parameter plane, it is unstable. For a weak Gamma distribution, the region of instability is constrained to a small region of the parameter plane with very small values of α and γ, and endemic equilibrium is stable in the rest of the parameter plane. In the case of a strong Gamma distribution, the situation is somewhat reminiscent of the case of discrete delay in that endemic equilibrium is unstable for smaller values of γ and stable for larger values of γ, but a major difference is that now the region of instability is quite small and constrained to very small values of γ. Finally, for the AF kernel, the picture is similar to that for the strong kernel in terms of an endemic steady state being unstable for smaller values of γ and stable for higher values of γ. The difference from the case of a strong Gamma distribution is that the value of γ, at which the transition from instability to stability takes place, decreases with α, as opposed to a strong Gamma distribution, where it was increasing with α.

FIG. 5.

Regions of stability and instability of the endemic steady state E 3 of model (12) depending on parameters α and γ. (a) Discrete delay kernel. (b) Weak Gamma distribution. (c) Strong Gamma distribution. (d) AF kernel. The color code denotes max[Re( λ)], with stable regions shown in blue-green and unstable shown in yellow-red. White indicates that the endemic steady state E 3 is biologically infeasible. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, τ = 100 days, k = 400.

FIG. 5.

Regions of stability and instability of the endemic steady state E 3 of model (12) depending on parameters α and γ. (a) Discrete delay kernel. (b) Weak Gamma distribution. (c) Strong Gamma distribution. (d) AF kernel. The color code denotes max[Re( λ)], with stable regions shown in blue-green and unstable shown in yellow-red. White indicates that the endemic steady state E 3 is biologically infeasible. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, τ = 100 days, k = 400.

Close modal

To better understand the effect of public heath information campaigns on (in)stability of an endemic steady state, in Fig. 6, we have plotted regions of stability of an endemic steady state for all four delay distributions and different values of γ. Compared to the situation shown in Fig. 2, which corresponds to the case γ = 0, as the impact of a public health information campaign becomes positive, we observe that for all delay distributions, regions of stability of endemic equilibrium increase with increasing values of γ. Counter-intuitively, this means that for the same imitation rates, an endemic steady state remains stable for higher values of mean time delay characterizing memory of past VSE. We also note that for the AF kernel, an increase in the size of the stability region for an endemic steady state is much more substantial than for the strong Gamma kernel, and it is not too dissimilar from a stability region observed for the case of a weak Gamma kernel. This suggests that more impactful public health information campaigns result in making the case of the AF kernel become more akin to the situation, where there is a discrete period of memory of past VSE.

FIG. 6.

Regions of stability and instability of the endemic steady state E 3 for different delay distribution kernels depending on the mean time delay τ and imitation rate k. Top row: discrete delay kernel, second row: weak Gamma distribution, third row: strong Gamma distribution, and fourth row: AF kernel with T 1 , 2 = ( 0.5 ± 0.4 ) τ. Left column: γ = 10 6, middle column: γ = 10 5, and right column: γ = 10 4. The color code denotes max[Re( λ)], with stable regions shown in blue-green and unstable shown in yellow-red. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, α = 0.002.

FIG. 6.

Regions of stability and instability of the endemic steady state E 3 for different delay distribution kernels depending on the mean time delay τ and imitation rate k. Top row: discrete delay kernel, second row: weak Gamma distribution, third row: strong Gamma distribution, and fourth row: AF kernel with T 1 , 2 = ( 0.5 ± 0.4 ) τ. Left column: γ = 10 6, middle column: γ = 10 5, and right column: γ = 10 4. The color code denotes max[Re( λ)], with stable regions shown in blue-green and unstable shown in yellow-red. Parameter values are μ = 3.9 × 10 5 days 1, v = 1 / 7 days 1, α = 0.002.

Close modal

In this paper, we have analyzed dynamics of pediatric vaccination when modeled as an imitation game, where a perceived risk of vaccine is determined as a delayed history of side effects, represented as an integral with some delay kernel. We have identified parameter regions of feasibility and stability of different steady states and obtained conditions for their stability. The results show that while the basic reproduction number that controls disease eradication only depends on parameters characterizing the disease, stability of the endemic steady state is determined not only by the mean time delays, but also by the shape of a specific delay distribution. Exponentially fading memory, which is often assumed in epidemic models, provides largest values of mean time delays, for which the endemic steady state remains stable for the same values of the imitation rate.

Since a major role in achieving high vaccination take-up is played not only by the spread of information through direct social contacts but also by targeted public health information campaigns, we have analyzed the effects such campaigns have on model dynamics. Increasing the impact of these campaigns leads to a reduction in the basic reproduction number, thus making it easier to achieve disease eradication, and this result is independent of the delay distribution characterizing memory of past VSE. Interestingly, and somewhat counter-intuitively, higher values of the intensity of public health information campaigns are associated with larger regions of stability of endemic equilibrium for all delay distributions. While in the absence of public health campaigns, dynamics of the model with an acquisition-fading kernel is more similar to that of the strong Gamma distribution, as the impact of public health information campaigns increases, the stability region of an endemic steady state for the AF kernel more closely resembles that for discrete delay, suggesting that the effect of memory of past VSE becomes sharper in that they are rather remembered for some fixed period of duration.

There are several directions, in which the work presented in this paper could be extended. One possibility is to include time-dependent public health information campaigns and determine optimal information strategies using techniques of control theory, as it was done in Buonomo et al.26,27 To achieve better biological realism, the model can be modified to include seasonal variation in the disease transmission rate, which is known to be a major factor, particularly, for childhood diseases.29–31 For many infectious diseases, there is a non-negligible period of time, where individuals are already exposed to disease and are incubating it, but are not yet infectious. In the model, this could be implemented by introducing an extra compartment of exposed individuals or by representing an infection process as some distribution, such as a Gamma distribution with a certain number of stages, which can vary from 3 for SARS32,33 or influenza34,35 to 20 for measles.36 Another interesting avenue for research would be to consider vaccination dynamics with delay-distributed memory of past VSE on social networks,37,38 which could provide a more realistic representation of how spread of information through social contacts affects vaccination decision-making. Finally, an interesting and important question in the context of modeling dynamics of infectious diseases is the role of stochastic effects. In the specific context of vaccination dynamics considered in this paper, three issues would be of particular interest: stochastic amplification, stochastic switching of strategies, and stochastic extinction. Stochastic amplification is a phenomenon, where sustained stochastic oscillations are observed in individual stochastic realizations around a deterministically stable endemic steady state.39,40 The model presented in this paper effectively describes mean-field dynamics for a very large population, and an alternative would be to allow payoff gains for vaccinator and unvaccinator strategies to rather determine probabilities of switching to an alternative strategy.41,42 Finally, stochastic extinction concerns an observation that if disease prevalence becomes very low during inter-epidemic periods [which is what can be seen, e.g., in Fig. 4(l)], due to stochasticity, disease can completely disappear from the population,43–45 and it would be insightful to explore the extent to which this can happen in the present model.

The authors are grateful to Jess Abbott, with whom we discussed some initial aspects of this work.

The authors have no conflicts to disclose.

Y. N. Kyrychko: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Visualization (equal); Writing – original draft (equal). K. B. Blyuss: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Visualization (equal); Writing – original draft (equal).

The data that support the findings of this study are available within the article.

1.
V.
Ramalingaswami
, “
Psychosocial effects of the 1994 plague outbreak in Surat, India
,”
Mil. Med.
166
,
29
30
(
2001
).
2.
J.
Savulescu
,
J.
Pugh
, and
D.
Wilkinson
, “
Balancing incentives and disincentives for vaccination in a pandemic
,”
Nat. Med.
27
,
1500
1503
(
2021
).
3.
S.
Loomba
,
A.
de Figueiredo
,
S. J.
Piatek
,
K.
de Graaf
, and
H. J.
Larson
, “
Measuring the impact of COVID-19 vaccine misinformation on vaccination intent in the UK and USA
,”
Nat. Hum. Behav.
5
,
337
348
(
2021
).
4.
M. S.
Islam
,
A.-H. M.
Kamal
,
A.
Kabir
,
D. L.
Southern
,
S. H.
Khan
,
S. M.
Hasan
,
T.
Sarkar
,
S.
Sharmin
,
S.
Das
,
T.
Roy
, and
M. G. D.
Harun
, “
COVID-19 vaccine rumors and conspiracy theories: The need for cognitive inoculation against misinformation to improve vaccine adherence
,”
PLoS One
16
,
e0251605
(
2021
).
5.
M. A.
Penţa
and
A.
Băban
, “
Mass media coverage of HPV vaccination in Romania: A content analysis
,”
Health Educ. Res.
29
,
977
992
(
2014
).
6.
R. D.
Todor
,
G.
Bratucu
,
M. A.
Moga
,
A. N.
Candrea
,
L. G.
Marceanu
, and
C. V.
Anastasiu
, “
Challenges in the prevention of cervical cancer in Romania
,”
Int. J. Environ. Res. Public Health
18
,
1721
(
2021
).
7.
P.
Fine
,
K.
Eames
, and
D. L.
Heymann
, “
`Herd immunity': A rough guide
,”
Clin. Infect. Dis.
52
,
911
916
(
2011
).
8.
C. J. E.
Metcalf
,
M.
Ferrari
,
A. L.
Graham
, and
B. T.
Grenfell
, “
Understanding herd immunity
,”
Trends Immunol.
36
,
753
755
(
2015
).
9.
B.
Ashby
and
A.
Best
, “
Herd immunity
,”
Curr. Biol.
31
,
R174
R177
(
2021
).
10.
R. M.
Anderson
and
R. M.
May
,
Infectious Diseases of Humans: Dynamics and Control
(
Oxford University Press
,
Oxford
,
1991
).
11.
F.
Guerra
,
S.
Bolotin
,
G.
Lim
,
J.
Heffernan
,
S. L.
Deeks
,
Y.
Li
, and
N. S.
Crowcroft
, “
The basic reproduction number ( R 0 ) of measles: A systematic review
,”
Lancet Infect. Dis.
17
,
e420
e428
(
2017
).
12.
P.
Plans-Rubió
, “
Low percentages of measles vaccination coverage with two doses of vaccine and low herd immunity levels explain measles incidence and persistence of measles in the European Union in 2017-2018
,”
Eur. J. Clin. Microbiol. Infect. Dis.
38
,
1719
1729
(
2019
).
13.
A.
d’Onofrio
,
P.
Manfredi
, and
E.
Salinelli
, “
Fatal SIR diseases and rational exemption to vaccination
,”
Math. Med. Biol.
25
,
337
357
(
2008
).
14.
C. T.
Bauch
, “
Imitation dynamics predict vaccinating behaviour
,”
Proc. R. Soc. B
272
,
1669
1675
(
2005
).
15.
D. A.
Salmon
,
S. P.
Teret
,
C. R.
MacIntyre
,
D.
Salisbury
,
M. A.
Burgess
, and
N. A.
Halsey
, “
Compulsory vaccination and conscientious or philosophical exemptions: Past, present and future
,”
Lancet
367
,
436
442
(
2006
).
16.
T. C.
Reluga
,
C. T.
Bauch
, and
A. P.
Galvani
, “
Evolving public perceptions and stability in vaccine uptake
,”
Math. Biosci.
204
,
185
198
(
2006
).
17.
A.
d’Onofrio
,
P.
Manfredi
, and
P.
Poletti
, “
The impact of vaccine side effects on the natural history of immunization programmes: An imitation-game approach
,”
J. Theor. Biol.
273
,
63
71
(
2011
).
18.
P. Y.
Geoffard
and
T.
Philipson
, “
Disease eradication: Private versus public vaccination
,”
Am. Econ. Rev.
87
,
222
230
(
1997
), see https://www.jstor.org/stable/2950864.
19.
A.
d’Onofrio
and
P.
Manfredi
, “
Vaccine demand driven by vaccine side effects: Dynamic implications for SIR diseases
,”
J. Theor. Biol.
264
,
237
252
(
2010
).
20.
P.
Manfredi
and
A.
d’Onofrio
,
Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases
(
Springer
,
New York
,
2013
).
21.
Z.
Wang
,
C. T.
Bauch
,
S.
Bhattacharyya
,
A.
d’Onofrio
,
P.
Manfredi
,
M.
Perc
,
N.
Perra
,
M.
Salathé
, and
D.
Zhao
, “
Statistical physics of vaccination
,”
Phys. Rep.
664
,
1
113
(
2016
).
22.
C. T.
Bauch
and
D. J.
Earn
, “
Vaccination and the theory of games
,”
Proc. Natl. Acad. Sci. U.S.A.
101
,
13391
13394
(
2004
).
23.
J.
Hofbauer
and
K.
Sigmund
,
Evolutionary Games and Population Dynamics
(
Cambridge University Press
,
Cambridge
,
1998
).
24.
A.
d’Onofrio
,
P.
Manfredi
, and
P.
Poletti
, “
The interplay of public intervention and private choices in determining the outcome of vaccination programmes
,”
PLoS One
7
,
e45653
(
2012
).
25.
B.
Buonomo
,
G.
Carbone
, and
A.
d’Onofrio
, “
Effect of seasonality on the dynamics of an imitation-based vaccination model with public health intervention
,”
Math. Biosci. Eng.
15
,
299
321
(
2018
).
26.
B.
Buonomo
,
P.
Manfredi
, and
A.
d’Onofrio
, “
Optimal time-profiles of public health intervention to shape voluntary vaccination for childhood diseases
,”
J. Math. Biol.
78
,
1089
1113
(
2019
).
27.
B.
Buonomo
,
R.
Della Marca
, and
A.
d’Onofrio
, “
Optimal public health intervention in a behavioural vaccination model: The interplay between seasonality, behaviour and latency period
,”
Math. Med. Biol.
36
,
297
324
(
2019
).
28.
A.
d’Onofrio
,
P.
Manfredi
, and
E.
Salinelli
, “Vaccinating behaviour and the dynamics of vaccine preventable diseases,” in Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases, edited by P. Manfredi and A. d’Onofrio (Springer, New York, 2013), pp. 267–287.
29.
S.
Altizer
,
A.
Dobson
,
P.
Hosseini
,
P.
Hudson
,
M.
Pascual
, and
P.
Rohani
, “
Seasonality and the dynamics of infectious diseases
,”
Ecol. Lett.
9
,
467
484
(
2006
).
30.
K. M.
Bakke
,
M. E.
Martinez-Bakker
,
B.
Helm
, and
T. J.
Stevenson
, “
Digital epidemiology reveals global childhood disease seasonality and the effects of immunization
,”
Proc. Natl. Acad. Sci. U.S.A.
113
,
6689
6694
(
2016
).
31.
B.
Buonomo
,
N.
Chitnis
, and
A.
d’Onofrio
, “
Seasonality in epidemic models: A literature review
,”
Ric. Mat.
67
,
7
25
(
2018
).
32.
C. T.
Bauch
,
J. O.
Lloyd-Smith
,
M. P.
Coffee
, and
A. P.
Galvani
, “
Dynamically modeling SARS and other newly emerging respiratory illnesses: Past, present, and future
,”
Epidemiology
16
,
791
801
(
2005
).
33.
S.
Riley
,
C.
Fraser
,
C. A.
Donnelly
,
A. C.
Ghani
,
L. J.
Abu-Raddad
,
A. J.
Hedley
,
G. M.
Leung
,
L. M.
Ho
,
T. H.
Lam
,
T. Q.
Thach
, and
P.
Chau
, “
Transmission dynamics of the etiological agent of SARS in Hong Kong: Impact of public health interventions
,”
Science
300
,
1961
1966
(
2003
).
34.
M. J.
Keeling
and
P.
Rohani
,
Modeling Infectious Diseases in Humans and Animals
(
Princeton University Press
,
Princeton, NJ
,
2008
).
35.
S.
Cauchemez
,
F.
Carrat
,
C.
Viboud
,
A. J.
Valleron
, and
P. Y.
Boelle
, “
A Bayesian MCMC approach to study transmission of influenza: Application to household longitudinal data
,”
Stat. Med.
23
,
3469
3487
(
2004
).
36.
R. E.
Hope-Simpson
, “
Infectiousness of communicable diseases in the household (measles, chickenpox, and mumps)
,”
Lancet
260
,
549
554
(
1952
).
37.
F.
Fu
,
D. I.
Rosenbloom
,
L.
Wang
, and
M. A.
Nowak
, “
Imitation dynamics of vaccination behaviour on social networks
,”
Proc. R. Soc. B
278
,
42
49
(
2011
).
38.
M. L.
Ndeffo Mbah
,
J.
Liu
,
C. T.
Bauch
,
Y. I.
Tekel
,
J.
Medlock
,
L. A.
Meyers
, and
A. P.
Galvani
, “
The impact of imitation on vaccination behavior in social contact networks
,”
PLoS Comput. Biol.
8
,
e1002469
(
2012
).
39.
D.
Alonso
,
A. J.
McKane
, and
M.
Pascual
, “
Stochastic amplification in epidemics
,”
J. R. Soc. Interface
4
,
575
582
(
2007
).
40.
R.
Kuske
,
L. F.
Gordillo
, and
P. E.
Greenwood
, “
Sustained oscillations via coherence resonance in SIR
,”
J. Theor. Biol.
245
,
459
469
(
2007
).
41.
T.
Maruta
, “
Binary games with state dependent stochastic choice
,”
J. Econ. Theory
103
,
351
376
(
2002
).
42.
I.
Soo Lim
and
P.
Wittek
, “
Satisfied-defect, unsatisfied-cooperate: An evolutionary dynamics of cooperation led by aspiration
,”
Phys. Rev. E
98
,
062113
(
2018
).
43.
Y.
Cai
,
Y.
Kang
,
M.
Banerjee
, and
W.
Wang
, “
A stochastic epidemic model incorporating media coverage
,”
Commun. Math. Sci.
14
,
893
910
(
2016
).
44.
B.
Cao
,
M.
Shan
,
Q.
Zhang
, and
W.
Wang
, “
A stochastic SIS epidemic model with vaccination
,”
Phys. A
486
,
127
143
(
2017
).
45.
Y.
Cai
,
J.
Jiao
,
Z.
Gui
,
Y.
Liu
, and
W.
Wang
, “
Environmental variability in a stochastic epidemic model
,”
Appl. Math. Comput.
329
,
210
226
(
2018
).
Published open access through an agreement with University of Sussex