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.

## I. INTRODUCTION

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\u22121/ 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 18^{10} (with some more recent estimates suggesting a range of between 6.51 and 45.4 in European countries^{11}), 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’Onofrio^{20} 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 Bauch^{14,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 transmission^{26} or seasonally forced transmission.^{27}

^{19,28}Another approach, which was briefly mentioned in d’Onofrio

*et al.*,

^{17}uses the imitation game model

^{14}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:

^{14,17,24}we will assume perceived risk of infection $ r i$ to be proportional to disease prevalence; i.e., $ r i=\omega I$. Perceived risk of vaccination $ r v$ is taken in the form that depends on past vaccine-related side effects,

^{17}) and $n=2$ (strong kernel) turns into $ g ( s ) = \sigma e \u2212 \sigma s$ and $ g ( s ) = \sigma 2 s e \u2212 \sigma s$, respectively, and an acquisition-fading (AF) kernel

^{28}$ g ( s ) = ( e \u2212 s / T 1 \u2212 e \u2212 s / T 2 ) / ( T 1 \u2212 T 2 )$. In the limit $ T 1\u21920$, $ T 2=T$ (or $ T 2\u21920$, $ T 1=T$), this kernel transforms into the weak kernel with $\sigma =1/T$, and in the limit $ T 1\u2192 T 2=T$, AF kernel becomes the strong kernel with the same $\sigma =1/T$. These delay kernels are illustrated in Fig. 1, with the exception of discrete delay kernel given by the Dirac $\delta $-function at $ \tau 0$. The mean time delay for each of these delay distributions can be found as

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.

## II. STEADY STATES AND THEIR STABILITY

^{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 $\beta >(\mu +v)$.

## III. NUMERICAL BIFURCATION ANALYSIS AND SIMULATIONS

To explore how stability of the endemic steady state $ E 3$ changes with parameters, we fix the values of $\mu =1/(70\u22c5365)$ days $ \u2212 1$, which corresponds to an average life expectancy of 70 years, and $v=1/7$ days $ \u2212 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 $\beta $ through $\beta = R 0(\mu +v)$. Figure 2 illustrates how stability of endemic steady state $ E 3$ varies with the mean time delay $\tau $ 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).

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 $\tau $, 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\xb10.4)\tau $. 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 $\tau $ 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 $\tau $, 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 $\tau $ such that for a higher value of $\tau $, 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.

With the AF kernel providing a family of kernels with the same mean delay $\tau = 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 $\tau /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=(\tau \xb1W)/2$. Stability boundaries of endemic equilibrium in the $\tau $– $W$ plane show that increasing the width $W$ of the AF distribution results in increasing the range of mean time delays $\tau $, for which endemic equilibrium is stable. Interestingly, increasing the rate $\alpha $, 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 $\tau $. 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.

## IV. INFLUENCE OF PUBLIC HEALTH INFORMATION CAMPAIGNS

*et al.*

^{24}and introduce an extra term in equation for $ p \u02d9$, namely,

*et al.*

^{25}to show that for $\gamma \u2265\alpha $, 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:

In Fig. 5, we have fixed the value of imitation rate $k$ and the mean time delay $\tau $ and plotted regions of stability and instability of endemic equilibrium $ E ~ 3$ of the modified model (12) depending on the strength $\gamma $ of the public health information campaign and the rate $\alpha $, 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 $\gamma $, 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 $\alpha $ and $\gamma $, 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 $\gamma $ and stable for larger values of $\gamma $, but a major difference is that now the region of instability is quite small and constrained to very small values of $\gamma $. 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 $\gamma $ and stable for higher values of $\gamma $. The difference from the case of a strong Gamma distribution is that the value of $\gamma $, at which the transition from instability to stability takes place, decreases with $\alpha $, as opposed to a strong Gamma distribution, where it was increasing with $\alpha $.

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 $\gamma $. Compared to the situation shown in Fig. 2, which corresponds to the case $\gamma =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 $\gamma $. 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.

## V. DISCUSSION

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 SARS^{32,33} or influenza^{34,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.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Infectious Diseases of Humans: Dynamics and Control*

*Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases*

*Evolutionary Games and Population Dynamics*

*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.

*Modeling Infectious Diseases in Humans and Animals*