Vaccination games and imitation dynamics with memory

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.¹ 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,⁵ resulted in a very low take-up of vaccination, leading to Romania having the highest rates of cervical cancer in Europe.⁶ 

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 18¹⁰ (with some more recent estimates suggesting a range of between 6.51 and 45.4 in European countries¹¹), 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,¹² 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.¹⁸ 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²⁰ provides an excellent overview of many techniques used in such modeling (see also a more recent review paper by Wang et al.²¹). 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.²³ d’Onofrio et al.²⁴ 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.²⁵ 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²⁶ or seasonally forced transmission.²⁷ 

The work by d’Onofrio et al.¹⁷ 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.

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

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.

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 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 $α$⁠.

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.

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.³⁶ 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.