Spatio-temporal games of voluntary vaccination in the absence of the infection: the interplay of local versus non-local information about vaccine adverse events

1.   Introduction

Mathematical epidemiology (ME) of infectious disease is one of the oldest and richest areas of mathematical biology ^[1,2], but it is also the area that probably had the largest impact on actual public health policies. Indeed, mathematical models of infections are nowadays routinely applied by international and national public health institutions: from the design of immunization programs up to the preparedness plans against possible future pandemics ^[3].

A limitation of classical ME is that its key assumptions are the legacy of classical Statistical mechanics. As a consequence, the social contacts between agents at risk of spreading or acquiring the infection are treated as 'encounters of particles' of a perfect gas. Pairwise, the contagion process is abstracted as it were a chemical reaction. As a consequence, transmission processes have been traditionally modelled by means of the mass action law.

Treating human beings as molecules implies that the spread of infectious diseases is totally unaffected by the agents' behavior, and vice-versa, simply because behavior is absent from the models of classical ME. This means that e.g., models used to evaluate the impact and cost-effectiveness of vaccination programs under voluntary immunization did not include neither the individuals risk perceptions about the disease nor those about vaccine-related adverse events (VAEs). Similarly, ME models treat the social contact process between individuals as a physical constant, implying that individuals continue contacting each other at the same rate independently of the magnitude they perceived of the risk of contracting the infection.

Nowadays we know well that the above assumption are quite coarse, possible useful to depict 'normal' situations but totally inadequate to describe scenarios such as vaccine scares, where the perceived risk of VAEs blows up for a while, or the course of an epidemics with high mortality ^[4,5].

The pioneering work that first included the human behavior in ME was ^[6], which extended the classical Kermack and McKendrick's ODE epidemic model to account for behavioral responses. Unfortunately, this paper remained relatively isolated until the need of embedding human behavior in ME models became increasingly urgent. This need arose due to the onset of a range of new phenomena such as the increasing mistrust and opposition towards vaccines ^[7,8,9,10]. In ^[4,5] it was argued that such phenomena are characteristic of the current landscape of infection and public health in modern industrialised countries. Notably, mistrust towards vaccination can be considered as part of the more global phenomenon known as post-trust society ^[11].

This led in the last two decades to the birth of a new branch of ME: the behavioral Epidemiology of infectious diseases (BEID) ^[4,5]. The main aim of BEID is to properly model the role of human behavior in the transmission and control of infectious diseases. This is done by integrating the classical ME tools with models and ideas from disciplines ranging from psychology to neural sciences, and from economics to sociology.

Most BEID models of vaccination behavior include the strategic behavior of agents, and therefore extensively use Game Theory. In particular, Bauch ^[12] pioneered the applications of evolutionary games to describe the dynamics of the vaccine propensity in a population by using an Imitation Game Dynamics (IGD) i.e., a model where the strategy perceived as better at a given time spreads in the population through imitation. In ^[12] the perceived risk of infection is taken as linearly increasing with the infection prevalence, whereas the perceived risk of VAE is assumed constant. Instead, in ^[13] the perceived risk of suffering a vaccine side effect is modelled as an increasing function of the information of the present and past incidence of VAEs.

As expected, in BEID a key role is played by the information that agents can access and use to evaluate risks, which modulates their vaccine–related decisions. In ^[14], the IGD vaccination model was further extended by including the effect of public interventions in communicating actual risks from vaccines and infection.

Most previously cited BEID models ^{[12,13,14,15]}, and other works in the same line, are based on ordinary differential equation (ODE) models, and disregard any type of structural heterogeneity. In particular, they are 'spatially homogeneous', which is a crude approximation from at least two standpoints: (ⅰ) the wide and complicate population mobility patterns; (ⅱ) the complicate role played by the spatial information network for behavior in relation to health.

In relation to this, a still unexplored area in the BEID literature regards the incorporation of behavioral hypotheses within classical PDE models of spatial dynamics of the reaction-diffusion type. In our opinion this represent a worthwhile effort for two main reasons. The first one is substantive and relates to the critical role played by information in behavioral epidemiology models. From this standpoint it is fundamental to subdivide the type of information that can be accessed by vaccination decision makers into a few sharply distinct types namely, 'local' vs 'global' vs 'non-local' based on simple hypotheses on the underlying spatial information kernels. The second reason relates to the robust analytic techniques that classical diffusion mathematics makes it available for the understanding of real world processes. Indeed, since the pioneering works of Kendall ^[16] and of Bailey ^[1,17] published in sixties, mathematical reaction-diffusion theory was extensively applied in ME models for the spread of human ^{[2,18,19,20,21,22]} and animal ^[23] infectious diseases. A paradigmatic example is the well–known pioneering paper by Noble ^[24] on European plague epidemics in the 14th century. Recently, Ducrot and Giletti ^[25] showed that, under Fickian diffusion hypothesis, that the Kermack– McKendrick epidemic model with non-diffusive susceptible population can have pulsating traveling wave solutions. Very recently Magal and coworkers ^[26] adopted Fickian diffusion with anisotropic diffusion coefficients to model the spread of infectious diseases, with focus on the impact of diffusion on the basic reproduction number. They adopted a similar approach ^[27] to model the spread of influenza in Puerto Rico, including also behavioral effects. They obtained a good match between their simulations and available spatiotemporal data. As far as vector-borne infectious diseases are concerned, recently Fitzgibbon et al. ^[28] proposed a model, where the diffusion of hosts was described by Fickian diffusion. Zhao and colleagues ^[29] again adopted Fickian diffusion to model the spread of a two–groups infectious diseases of SIR type, focusing in particular on the onset of traveling waves.

As stressed in the review paper ^[30], in its pioneering period, an important chapter of the spatial modelling of the spread of infectious disease was represented by the modelling of the contagion as a non–local process, as in ^{[1,2,16,17,18,19,20]}.

In the light of the above background, our aim in this article is to investigate the interplay between vaccination dynamics, human decisions, spatial mobility and information. To do so we extended to a continuous spatially structured setting the behavioral epidemiology models for the impact of vaccine adverse events (VAEs) on vaccination choices that were first introduced in ^[13,14] by a classical economy–oriented game theory approach and reformulated as a process of double contagion of ideas in ^[5,15,31]. We note here that the first phenomenological approach is partially remindful of the phenomenological theory of innovation diffusion (TID) by Mahajan ^[32,33,34] (see also ^[36]). However, taking into the account the second approach, which is more mechanistic, important differences emerge with TID, which are stressed in Section 12.

In particular, we will consider voluntary immunization decisions and focus on the simple but important case where the disease is absent in the population under study. This case is not special at all. On the contrary, it is of central relevance since, e.g., it represents the case of a population where a previously endemic vaccine preventable infectious disease of childhood has been successfully eliminated, so that prevalence is equal to zero, but there is the need to maintain a high-coverage immunization policy in the post-elimination period to prevent the risk of infection re-emergence. A major instance is that of poliomyelitis in industrialized countries. In such a context, where the absence of the infection will remarkably reduce the incentive to immunize thereby weakening the probability of switching from the 'non-vaccinator' to the 'vaccinator' strategy, the dynamics of VAEs will arguably become the key determinant of vaccination decision and collective coverage. To investigate the interplay with spatial mobility the resulting vaccination dynamics is set into the simplest possible framework for spatial mobility, namely classical diffusion based on Fick's law.

The consideration of the spatial effects i.e., the spatial distribution of VAEs, requires to carefully take into account the information on VAEs that is handled by parents of children eligible for vaccination while forming their perceptions of risk, that they will subsequently use to take their immunization decisions.

In relation to this 'information issue' we considered three main scenarios. In the first scenario, the information that individuals access and use uniquely concerns VAEs occurred locally. In the second case, the information used is both the local one and a global, nation-wide, average. Roughly speaking, this scenario corresponds to the case where agents take their decisions also based on national media such as national news or the internet. The third scenario is an intermediate one and aims at taking into the account both the known phenomenon of information attenuation ^[37], and the reasonable assumption that agents give less weight to events occurred at far distant locations. We represented the latter scenario, where the information used is not purely local but it is also not global, by resorting to a range of different spatial information kernels.

Of the resulting wide collection of problems we have investigated (by also distinguishing the nature of space as bounded vs non-bounded) (ⅰ) the presence and stability of space homogeneous equilibria, focusing on the nontrivial behavioral equilibrium; (ⅱ) conditions for bifurcations; (ⅲ) existence of classical and generalized traveling waves; (ⅳ) effects of awareness campaigns enacted by the Public Health System to effectively sustain vaccine uptake, as first proposed by ^[14].

2.   The spatio-temporal imitation game for vaccination

To model the impact of human responses to VAEs on the overall dynamics of vaccine uptake in the new-born babies, we follow ^[12,13,14] and assume that vaccination is a binary decision problem, i.e., one for which only two, mutually exclusive, strategies are played by parents: "'vaccinator"' (strategy $1$), and "'non-vaccinator"' (strategy $2$). Let us consequently denote by $P(x, t)$ and $A(x, t)$ the fractions of vaccination decision makers (e.g., parents of children eligible for vaccination) at location $x$ and time $t$ that follow, respectively, strategy $1$ and $2$. The previously cited models of IGD ^[12,13,14] all relied on the concept of payoff. However, many evolutionary game models belong to the class of urn models ^[38], a remarkably wide family that includes both chemical kinetics models as well as classical ME models. Thus, we follow here ^[5,15,31], where the vaccination IGD was derived as a model of 'double contagion' of ideas between the two involved groups namely, the group playing the 'vaccinator' strategy and the one playing the 'non-vaccinator' strategy).

Thus, by extending to the present spatio-temporal setting the IGD in ^[5,31], we obtain the following family of models:

where $t>0$, $x\in\Omega$, with $\Omega$ a bounded subset of $\mathbb{R}^n$, $n = 1, 2$. If $n = 1$, then we assume that $\Omega = [-L, L]$, $L>0$, while, if $n = 2$, then we assume that its boundary is sufficiently smooth.

In particular, $\vartheta(M_{i})$ and $\alpha(M_{se})$ represent the strategy-specific 'transmission rates' following social contacts with individuals playing the other strategy. They are assumed to be non-decreasing functions of the variables $M_{i}$ and $M_{se}$. The latter are information indices summarizing the available information about the (current or past) state of the infection ($M_{i}$) and of VAEs ($M_{se}$), respectively, that are used by parents to formulate evaluations of related risks.

Assuming a stationary population and normalizing w.r.t. its steady state, since $A = 1-P$, system (2.1) yields:

As for the key quantities $M_i$ and $M_{se}$, given our focus on situations where the infection has been eliminated, we assume that the perceived risk related to infection is constant:

The previous formulation mimics the situation where the infection is absent, so the corresponding perceived risk is prevalence–independent. We assume that $\vartheta_{0}>0$ to reflect a non–null perceived risk of infection even in the absence of infection. This can be justified by the continued activity of an active public health system that aims to keep a high degree of population awareness on the risk of reintroduction which is perceived as homogeneous throughout the entire space. Moreover, we assumed that the perceived risk of VAEs depends on information on vaccine side effects occurring both locally but also non-locally, as follows:

where

where function $\phi(x)$ is non-negative with

and even:

To sum up, the perceived risk of VAEs $\alpha(M_{se})$, has three components: (ⅰ) a baseline value $\alpha_0$ possibly constant over space, mirroring a true underlying risk of VAEs, (ⅱ) a strictly local component $\alpha_1 P$, possibly reflecting local deviations in vaccine coverage and ensuing deviations in the local number of VAEs, and (ⅲ) a non-local component $\alpha_2 J(P)$ tuning risk perceptions arising from differences in VAEs at different geographic sites, according to a suitable spatial information kernel. In particular, the functional $J(P)$ describes two mathematically equivalent scenarios, where individuals located at position $x$ can: ⅰ) receive information concerning the value of $P$ at the space point $y$ with a weight/attenuation $\phi(x-y)$; ⅱ) receive full spatial information but, in taking their decisions, they assign a weight $\phi(x-y)$ to $P(y, t)$.

Setting $\vartheta_{*} = \vartheta_{0}-\alpha_0>0$, brings to the following model

with $x\in\Omega$ and $t>0$.

Let $l$ be a characteristic length of the domain $\Omega$; to mathematically simplify the analysis, the model is reformulated in terms of the following variables:

where we assumed $1/\vartheta_{*}$ is time unit. The set derived from $\Omega$ by adimensionalization will be denoted as $\Omega^{*}$.

Thus, the model (2.8) becomes (for the sake of simplicity we omit the stars):

with $x\in\Omega$ and $t>0$.

Moreover, by imposing Neumann boundary conditions:

where $n$ is the outer normal vector with respect to $\partial\Omega$.

In the case in which the characteristic size of $\Omega$ is very large (for example in one dimension, $L\gg 1$), then we can redefine $\Omega$ as $\Omega = (-L, L)$ and (as if we operated the limit $L\rightarrow +\infty$) assume that the domain is unbounded, i.e., $\Omega = \mathbb{R}^n$, $n = 1, 2$.

3.   The Fisher-Kolmogorov model as a particular case of (2.10)

In the special case where people do not take into account the spatial variation in VAEs, i.e., $\alpha(M_{se}(P)) = \alpha_0$, model (2.10) reduces to the well-known Fisher-Kolmogorov equation

introduced in ^[39,41]. The FK equation (3.1) in our case describes the spread (via traveling waves ^[22]) of the idealized case $P = 1$, where all people is in favor of the vaccination and the 'recession' of the case $P = 0$, where no subjects are in favor of vaccination.

4.   Modelling risk perceptions: global vs local vs non-local information

Arguably, the behavior of (2.10) critically depends on the specific functional form of the spatial information kernel $\phi(x)$. We list a number of relevant forms.

The simplest case is when the available information is purely local, that is

where $\delta(x)$ is the Dirac function. Then, Eq (2.10) reads as:

where, with slight abuse of notation, we set

As argued in the introduction, by far the most important specific form for $\phi(x)$ is the constant one, namely

where $\mu(\Omega)$ is the measure of $\Omega$; thus, the non–local term (2.5) yields:

This specific form models the noteworthy scenario where people are exposed to a global (average) information e.g., by the media (including the internet) and the Public Health system, about the current state of vaccine-side effects coming from the whole domain $\Omega$. The resulting model is:

Another important scenario is when the information about VAEs comes from the whole $\Omega$ but it is either attenuated with the distance (this was possibly true for diseases in historical epochs) or it is the individual receiving the information that pays lesser and lesser weight depending on the distance.

Two reasonable functional forms that might well represent this scenario are:

ⅰ) Gaussian decay

where $a>0$ and $1/\sqrt{a}$ is the characteristic attenuation length. In $\mathbb{R}^n$ ($n = 1, 2$) the normalization constant $C_g$ is:

In this case, model (2.10) becomes:

ⅱ) Exponential decay,

where $a>0$ and in $R^n$ ($n = 1, 2$) the constant $C_{e}$ reads as follows

In this case, model (2.10) becomes:

Another important class of spatial kernels is represented by non–local kernels with bounded support, i.e., that are null beyond a given distance. An example is:

where $Hev(.)$ is the Heaviside function and

In this case, model (2.10) becomes:

Another example in $\mathbb{R}$ is

and in $\mathbb{R}^2$

In this case, model (2.10) becomes:

As it is to be expected, the Fourier transform will be important in the mathematical analysis of the problem in study.

We use the following definition of Fourier transform for a generic function $f(x):\mathbb{R}^{n}\longrightarrow\mathbb{R}$ ($n = 1, 2$):

where $\xi$ is the Fourier's vector, $\xi x = \xi_1 x$ in $\mathbb{R}$ or $\xi x = \xi_1 x_1+\xi_2 x_2$ in $\mathbb{R}^2$. Here and further down $\xi^2 = \xi^2_1$ in $\mathbb{R}$ or $\xi^2 = \xi^2_1+\xi^2_2$ in $\mathbb{R}^2$.

We have that:

5.   Properties of the spatially homogeneous model

We briefly summarize the results in ^[13], concerning the behavior of the spatially homogenous model corresponding to (4.1). This is described by the ordinary differential equation (ODE)

A first equilibrium $P_0 = 0$ corresponds to the no vaccination scenario where no parent is favorable to immunization. This non-vaccinator equilibrium, NVE is always unstable. If $1 < \alpha$, then there exists a non-trivial behavioral equilibrium

that is globally attractive in (0, 1). We term $P_2$ a behavioral equilibrium (BE) because its level is tuned by the balance between the parents' perceptions of risks based on the handling of available information. Finally, there is a pure vaccinator equilibrium, PVE $P_1 = 1$ ^[12,13], corresponding to the unrealistic scenario where all parents are favorable to immunization. The PVE is unstable if $1 < \alpha$, and globally attractive if $0 \le \alpha \le 1$.

Arguably, in the case of an infection that has been eliminated one expects that at least in an initial phase it holds $P_2 > p_{C}$ where $p_{C} = 1-1/R_0$ is the critical vaccination threshold allowing infection elimination, and $R_0$ is the basic reproduction number of infection. However, this does not need to be true at subsequent times, and this is why it is important for the public health system to monitor the evolution of VAEs and related information. Finally, if $\alpha = 0$ then the model reduces to the well-known logistic model (see also next section).

6.   Stability analysis of homogeneous steady state

The equilibria of the ODE model (5.1) described in Section 5 are also homogenous equilibria of the spatially structured model (4.1). In this section we will investigate their stability, focusing primarily on the behavioral equilibrium $P_2$ defined in (5.2), since the other two equilibria ($P_0 = 0$, i.e., no vaccinators, and $P_1 = 1$, i.e., all vaccinators) do not correspond to realistic situations.

We will explore multiple relevant cases depending on: ⅰ) the nature of the information used to make the vaccination decision: local, global or non–local i.e., on the structure of the spatial information kernel; ⅱ) the nature of $\Omega$: bounded vs non-bounded.

6.1.   Local information

Under purely local information i.e., Eq (4.1) with boundary condition (2.11), the structure of the local stability of the homogeneous equilibria is similar to the one characteristic of the ODE model (5.1). We proved the following:

Theorem 1. The non-vaccinator equilibrium $P_0$ is always unstable.

If $\alpha>1$, the behavioral equilibrium $P_{2}$ is locally asymptotically stable (LAS) and the PVE $P_{1}$ is unstable. Moreover, if $\Omega$ is bounded then $P_2$ is also globally asymptotically stable (GAS).

If $0\le \alpha \le 1$, $P_{2}$ is not an admissible equilibrium and the PVE $P_{1}$ is LAS. Moreover, if $\Omega$is bounded then $P_1$ is also globally asymptotically stable (GAS).

We report the proof for the non-trivial parts. If $\Omega$ is bounded, from the linearized equation around the behavioral equilibrium $P_2 = 1/\alpha$,

it follows that the the m–th eigenvalue is given by:

where $\nu_m = -\xi_m^2$ is the m–th eigenvalue of the heat equation with Neumann conditions (2.11) in $\Omega$. For example, if $\Omega = [-L, L]$, then $\xi_m = m\pi/L$, $m\in\mathbb{N}$.

Similarly, if $\Omega$ is unbounded, then by applying the Fourier's Transform (4.4), we obtain the following expression for the spectrum:

Formulas (6.2) and (6.3) will be useful to assess the impact of the local vs non-local vs global structure of information on our system.

The global stability of $P_2$ and $P_1$, respectively, can both be demonstrated by adopting as Liapunov functional the free energy:

where $U(P)$ is the 'potential' (associated to the 'force' $P(1-P)(1- \alpha P)$), given by:

from which the claim easily follows.

6.2.   Global information

Under global information i.e., when people are exposed to the average information on VAEs coming from the whole domain (note this scenario only concerns the case where $\Omega$ is bounded):

the model (2.10) becomes:

with Neumann conditions (2.11).

The linearized equation at $P_2 = 1/\alpha$ reads as follows

Setting $w(x) = \cos(\xi_m x)$ as eigenfunction, it is straightforward to show that the m–th corresponding eigenvalue is given by

Thus $P_2$ is LAS.

Interestingly, if $m = 0$, referring to (6.2) and (6.7), then $\lambda_0 = \widetilde{\lambda}_0$. On the contrary, if $m\geq 1$, comparing eigenvalues in (6.2) and (6.7), we can deduce that $\widetilde{\lambda}_m>\lambda_m$: therefore eigenvalues in (6.7) are 'less negative' than their counterparts for the local information model in (6.2). This indicates a slower convergence to $P_2$ in Eq (6.5) compared to Eq (4.1) for all modes $m \ge 1$.

As for the NVE $P_0 = 0$, the eigenvalues obey

i.e., at least the mode $0$ is unstable. Moreover, the linearized equation of (6.5) around the PVE $P_1 = 1$ gives the eigenvalues:

Thus $P_1$ is unstable, unless the unlikely event $0\le \alpha \le 1$ holds, implyig the LAS of $P_1$.

6.3.   Non–local information

At variance with the global information scenario, under non-local information it is necessary to distinguish the case where $\Omega$ is unbounded from the one where it is bounded.

6.3.1.   Stability of homogeneous equilibria: $\Omega$ unbounded

As noted above, assuming that $L\gg 1$, we can study the following equation deduced from (2.10)

with $x\in\Omega = \mathbb{R}^n$, $n = 1, 2$, $t>0$ and Neumann condition (2.11) on the boundary (i.e., at infinity).

The associated linearized equation at a generic homogeneous steady state $P^{*}$ is:

The associated eigenvalue problem has the form:

Let us now focus on the case of primary interest here namely when $P^{*}$ is given by the behavioral equilibrium i.e., $P_2$ defined in (5.2). This yields:

We consider the linearized operator $\mathcal{L}$ on $L^{2}(\mathbb{R}^n)$ (n = 1, 2), such that

The spectrum of this operator is denoted by $\sigma(\mathcal{L})$. The eigenvalue problem

is called spectrally stable if $\sigma(\mathcal{L})\subset[0, +\infty)$ and spectrally unstable if there exists $\lambda <0\ |\ \lambda\in\sigma(\mathcal{L})$. Moreover, $\sigma(\mathcal{L})$ is spectrally stable iff the operator $\mathcal{L}$ is positive, i.e., $(\mathcal{L}W, W)>0$ for any $W$. Let us introduce the function

where $\tilde{\phi}(\xi)$ is the Fourier transform (defined in formula (4.4)) of $\phi(x)$, which is real since we assumed symmetric $\phi(x)$.

Noting that

we can state the following:

ⅰ) function $\Phi(\xi)$ evaluated at $\xi = 0$ is always positive, i.e.,

ⅱ) for any given continuous and bounded function $\tilde{\phi}$, $\Phi(\xi)$ becomes strictly positive if $D$ is sufficiently large. Thus, for sufficiently large values of the diffusion coefficient the behavioral equilibrium $P_2$ will be stable irrespective of the shape of the spatial information kernel.

On the contrary, if $\phi(x)$ has bounded support, then its Fourier transform can take negative values that can destabilize $P_2$.

We now consider the special cases of the spatial information kernel $\phi(x)$ considered in Section 4.

If $\phi(x) = \phi_1(x) = C_g e^{-ax^{2}}$ then function $\Phi$ takes the form:

As the latter expression is strictly positive, the behavioral equilibrium $P_{2}$ is always locally stable.

Similarly, under the exponential tent kernel $\phi(x) = \phi_2(x) = C_e e^{-a|x|}$, $x\in \mathbb{R}^n$ ($n = 1$ or $2$):

Therefore, also in this case function $\Phi(\xi)$ is strictly positive and $P_{2}$ is locally asymptotically stable.

These two examples suggest that kernels that are non-null over the whole space, and have a realistic shape (i.e., are decreasing in the distance from the current site), always promote the local stability of the behavioral equilibrium $P_2$.

To analyse the effects of the two proposed kernels with bounded support, we will only consider, for sake of notation simplicity, the case $\Omega = \mathbb{R}$.

Under the kernel $\phi(x) = \phi_3(x) = C_N (N^2-x^2) Hev(N-|x|)$, with $N>0$, function $\Phi$ is as follows:

Depending on $N$ and on the other parameters, $\Phi(\xi)$ can be strictly positive or can change sign. As noted before, $\Phi(0) = 1-1/\alpha >0$. Also in this case, if $\xi$ is large enough, the function $\Phi(\xi)$ is positive and the equilibrium $P_{2}$ is locally stable. When the diffusion coefficient is large enough, the local stability in ensured. The instability can appear if $D$ is small enough and if the ratio $\rho = \alpha_2/\alpha_1$ is large enough.

This type of behavior is illustrated in the Figure 1 in the particular case in which $\alpha_1 = 0$.

We finally consider $\phi(x) = \phi_4(x) = 1/(2h) Hev(h-|x|)$, with $h>0$. The function $\Phi(\xi)$ reads as follows:

Depending on $h$ and the other parameters, this function can be strictly positive or can change sign. Moreover, $\Phi(\xi)$ is positive for $|\xi|$ sufficiently large. If the diffusion coefficient is large, then $\Phi(\xi)$ is positive for all $\xi$, and the behavioral equilibrium $P_{2}$ is locally stable. The instability can appear if $D$ is small enough and $\rho = \alpha_2/\alpha_1$ is large enough.

A possible form of function $\Phi(\xi)$ is given in Figure 2, in the special case in which $\alpha_1 = 0$.

So we may say that in the examined examples, the stability of $P_2$ depends on whether or not the support of $\phi(x)$ is bounded. If it is bounded, then $P_2$ can become unstable depending on the parameters characterizing the kernel $\phi(x)$, i.e., on the features of the kernel. In this case some spatially heterogeneous solutions can appear.

Proceeding as above, it is an easy matter to show that replacing $P^{*} = P_0 = 0$ (the 'no-vaccinators' scenario) in eigenvalue problem (6.10), we obtain instability (since $\Phi(\xi) = D\xi^{2}-1$); moreover if $\alpha> 1$ (that correspond to biological existence of $P_2$), also $P_1 = 1$ is unstable:

In the unlikely event that $\alpha\leq 1$, instead, the all-vaccinator homogeneous equilibrium $P_1 = 1$ is LAS.

In Figure 3 an illustration of stability analysis of equilibrium $P_2$ is given for Eq (6.8). The Eq (2.10) is simulated with the step function $\phi_4(x) = 1/(2h)Hev(h-|x|)$ as kernel. Choosing a small perturbation as initial data, in Figure 3(a), the solution evolves in time approaching to the stable equilibrium $P_2$; while in Figure 3(b) the loss of stability of the constant equilibrium point $P_2$ and the convergence to a periodic spatial structure is shown. The numerical results confirm the analytical ones: in the case where $D$ is sufficiently large and the step $h$ small, $P_2$ is LAS, if $D$ is small and $h$ increases, we see convergence to non-constant structures.

6.3.2.   Stability of uniform equilibria: $\Omega$ bounded

Assuming now that $\Omega$ is bounded, we obtain the following equation from (2.10)

with $t>0$, $x\in\Omega$ and Neumann conditions (2.11).

For the sake of notation simplicity we will consider here $\Omega = [-L, L]$. The eigenvalues problem obtained as a result of the linearization about a generic steady state $P^*$ is given by the equality

Since for both the NVE and the PVE equilibria $P^* = P_{0} = 0$ and $P^* = P_{1} = 1$ the non–local term disappears, the resulting local stability analysis is straightforward and it is omitted. Therefore, we will only focus on the behavioral equilibrium $P_2 = 1/\alpha$.

Here, we make some assumptions. We will suppose that the function $\phi(x)$ is defined on the whole axis and that it is periodic with the period $2L$:

Taking $P^{*} = P_{2}$, brings to the eigenvalue problem

in the form

Note that the boundary conditions (2.11) are satisfied. Taking into account that function $\phi$ is periodic and even, we obtain

Thus, we have that the m–th eigenvalue reads as follows:

where

We note that $\lambda_{m}^{*}$ is negative for $m = 0$ and for $m$ sufficiently large, while it can be positive for some intermediate values of $m$ depending on the function $\phi$ and on the values of parameters, similarly to what we observed in the continuous spectrum analysis.

Let us now consider $\phi(x) = \phi_3(x) = (1/(2h))Hev(h-|x|)$, $0<h<L$. Then

Consider the functions

If $\Phi_{m}(h)$ is a positive function, then the corresponding eigenvalue $\lambda_{m}^{*}$ is also positive. Let us find the conditions on parameters when the maximal eigenvalue is zero (stability boundary). From the conditions

where prime denotes the derivative with respect to $h$, we obtain

where $\nu = \xi_m h$. The first relation in (6.25) allows us to find $\nu$, and the second relation determines the stability boundary.

7.   Bifurcation analysis for nonlocal information in case of bounded $\Omega$

In the case of nonlocal information and bounded one-dimensional $\Omega$, it is possible to carry out the bifurcation analysis (a.k.a. weakly nonlinear analysis) for Eq (6.19) under the assumption that $D$ is the bifurcation parameter ^[40].

The stationary non-homogenous solutions of the model solve

with $P'(-L) = P'(L) = 0$ and assuming that conditions (6.21) for the kernel $\phi(x)$ hold.

The equilibrium state $P(x) = P_2 = 1/\alpha$ is a solution of (7.1), independently from the values assumed by $D$.

If $D$ crosses a bifurcation value $D_0$, a simple real eigenvalue of the linearized problem crosses zero. In order to study this bifurcation, we look for solutions of (7.1) in the form of the expansion

where $\varepsilon$ is a small parameter. We set

Substituting these expansions into Eq (7.1) and equating the terms in $\varepsilon^{1}$, we get

with $p_{1}(-L) = p_{1}(L) = 0$. This problem coincides with eigenvalue problem (6.20) with $P^{*} = P_{2}$ and $\lambda = 0$. Hence the value $D_{0}$ should be chosen in such a way that this eigenvalue problem has a zero eigenvalue, $p_{1}(x) = \cos\left(\xi_m x\right)$ is the corresponding eigenfunction, with $m$ an integer, $m\neq 0$.

Next, we equate the terms with $\varepsilon^{2}$:

where $p_{2}(-L) = p_{2}(L) = 0$ and

In order to obtain solvability conditions for problem (7.3), let us note that problem (7.2) is self-adjoint since the kernel $\phi$ is an even function. Indeed, it can be directly verified that

where $\mathcal{L}$ is the operator which corresponds to the left-hand side of (7.2) and which acts on $C^{2}$ functions satisfying the boundary conditions. Hence problem (7.3) is solvable if and only if

Therefore

and

where

The terms with $\varepsilon^{3}$ give the problem

where

From the solvability condition

we obtain

If $D_{2}\neq 0$, then from the expansion for $D$ we obtain

and up to the second-order terms,

The bifurcation is supercritical ($D>D_{0}$) for $D_{2}>0$ and $P_2>\frac{1}{2}$. The bifurcation is subcritical ($D<D_{0}$) if the ratio $(2P_2-1)/D_2$ is negative.

8.   Traveling waves and generalized traveling waves

Preliminarily, we recall that in Section 2 we pointed out that in case of constant rate of transfer from the strategy 'vaccine' to the strategy 'no–vaccine' ($\alpha(M_{se}) = \alpha_0$) our model is equivalent to the most prototypical equation generating traveling waves, the Fisher-Kolmogorov equation (3.1). However, this case can more be considered as pathological than a realistic scenario.

In this section, we look in non-pathological cases (from the epidemiological viewpoint) for the onset of possible traveling waves ^[22,43], i.e., heteroclinic connections between two equilibria.

In this section, we restrict the analysis to the case in which the domain $\Omega$ is one-dimensional and unbounded, i.e., $\Omega = \mathbb{R}$.

8.1.   Local information

Wave solutions for Eq (4.1)

are function of the form $P(x-ct) = p(z)$ bounded on the whole axis and twice continuously differentiable. The constant $c$ is a speed of the wave, $z = x-ct$ is the moving coordinate frame and

with $P_{-}$ unstable and $P_{+}$ stable equilibrium point. After substitution of variable $z$ in Eq (4.1), we obtain the following ordinary differential equation:

It is known that, a necessary condition for the existence of the solution is that the $P_{-}$ must have an unstable (departing) manifold and $P_{+}$ a stable (incoming) manifold. After stability analysis, and following ^[44], we obtain the following result:

Theorem 2. If $\alpha>1$, there exists a constant $c^{(1)} = 2\sqrt{D}$ such that $\forall\ c\in[c^{(1)}, +\infty[$, there exists a traveling wave solution of velocity $c$ connecting equilibrium $P_{2} = 1/\alpha$ and equilibrium $P_{0} = 0$, i.e., a function $P(x-ct)$, solution of Eq (8.2) on the real line $]-\infty, +\infty[$ and satisfying Eq (8.1) with $P_{-} = P_{2}$ and $P_{+} = P_{0}$. This solution is monotone decreasing and the derivatives $P''(z)$ and $P'(z)$ tend to zero as $x\rightarrow\pm\infty$. Moreover, a traveling wave solution of velocity $c$ connecting equilibria $P_{1}$ and $P_{0}$ cannot exist.

We recall that for bounded $\Omega$ we had shown that the spatially homogenous equilibrium $P_2$ is GAS. In the above theorem it is shown, roughly speaking, that the TW are such that the scenario where the value $P = P_2$ spreads until it invades all the space.

The case $0<\alpha<1$ has here some further mathematical interest, since the equation can be read as a modification of the Fisher-Kolmogorov equation.

Theorem 3. If $0<\alpha<1$, there exists a constant $c^{(1)}$ such that $\forall\ c\in[c^{(1)}, +\infty[$, there exists a traveling wave solution of velocity $c$ connecting equilibrium $P_{1} = 1$ and equilibrium $P_{0} = 0$, i.e., a function $P(x-ct)$, solution of Eq (8.2) on the real line $]-\infty, +\infty[$ and satisfying Eq (8.1) with $P_{-} = P_{1}$ and $P_{+} = P_{0}$. This solution is monotone decreasing and the derivatives $P''(z)$ and $P'(z)$ tend to zero as $x\rightarrow\pm\infty$.

Traveling wave solutions for the Eq (4.1) are shown in Figure 4 in the case in which $\alpha>1$. They connect the constant state $P_2$ at $-\infty$ and the constant state $P_0$ at $+\infty$, with the minimal speed $c = c^{(1)} = 2\sqrt{D}$.

8.2.   Non–local information

In previous sections, we proved that, in the case of non–local information, the spatially homogeneous equilibrium solution $P_{2} = 1/\alpha$ can lose its stability. Then some spatial structures can bifurcate from it. Therefore, instead of traveling waves connecting $P_{0}$ and $P_{2}$, we can expect the existence of some other solutions connecting $P_{0}$ at $+\infty$ with some structures at $-\infty$. Such solutions are called generalized traveling waves (GTW) and were first introduced in ^[45] for reaction-diffusion systems. Such solutions can be characterized by two main properties:

(1) They exist for all $t\in\mathbb{R}$. Moreover, under some conditions, such solution can be unique and stable.

(2) They are propagating solutions, which can be explained as follows: let $q$ be a constant, $P_{+}<q<P_{-}$. For each $t$ fixed consider the equation $P(x, t) = q$ with respect to $x$. Denote by $m^{+}_{q}(t)$ its maximal solution (if it exists) and by $m^{-}_{q}(t)$ its minimal solution. If $\frac{m^{\pm}_{q}(t)}{t}\rightarrow c$ as $t\rightarrow\infty$, then we say that this solution propagates with the speed $c$. Thus, GTW are global propagating solutions.

Thus, we want to study GTW solution for the equation

where $\Omega$ is unbounded.

Consider the Cauchy problem for the equation

with $c\geq 2\sqrt{D}$ a given constant. Assume that the initial condition $P(x, 0) = P^{0}(x)$ is non-negative and less than 1. Then the solution $P(x, t)$ exists and is also non-negative and less than 1 for all $t\geq 0$. Consider

Hence

Equation (8.4) becomes

Consider also the equation

The following result holds:

Lemma 4. If $c\geq 2\sqrt{D}$, there exists a stationary solution $v_{c}(x)$ of (8.6) such that

Proof. Denoting $z = v-P$ and taking the difference of Eqs (8.6) and (8.5), we obtain

Since the last term in the right-hand side of Eq (8.8) is non-negative, then from the inequality $z(x, 0)>0$ for all $x\in\mathbb{R}$, it follows that $z(x, t)>0$ for all $t>0$ and $x\in\mathbb{R}$. Hence,

Stationary solutions of Eq (8.6) correspond to traveling waves for the KPP-equation ^[39,41]. If $c\geq 2\sqrt{D}$, then they are monotone in space and stable, while for $0<c<2\sqrt{D}$ they are non-monotone and unstable ^[41,44]. These waves have limits at infinity: $v(-\infty) = 1$, $v(+\infty) = 0$.

Denoting by $v_{c}(x)$ a stationary solution of Eq (8.6) for $c\geq 2\sqrt{D}$, the claim is proved.

Thus, we obtain an estimate from above of the solution of the Cauchy problem associated to (8.5). We now estimate it from below. From the inequality (8.9),

We assume $w<1$: consider the equation

Denoting $s = P-w$ and taking the difference between Eqs (8.5) and (8.11), we obtain

Since the last two terms in the right-hand side of this equation are non-negative, then

The following result holds:

Lemma 5. If $c\geq 2\sqrt{D}$, there exists a stationary solution $w_{c}(x)$ of (8.11) such that $w_{c}(x_{0}) = 0$ for some $x_{0}$, $w_{c}(x)>0$ for $x>x_{0}$, $w_{c}(x)<0$ for $x<x_{0}$ and $w_{c}(x)\sim v_{c}(x)$, when $x\rightarrow\infty$

Proof. We look for stationary solution of Eq (8.11), i.e.,

Since $v_{c}$ and $K(x)$ tend to zero when $x\rightarrow\infty$, then $1-\alpha_{1}v_{c}-\alpha_{2}K(x)$ is close to 1 in some right half-axis. Then, for $c\geq 2\sqrt{D}$, solutions of (8.14) are close to functions $v_{c}(x)$ when $x\rightarrow\infty$.

We prove the existence of a solution which has a zero. Denote $g(x) = 1-\alpha_{1}v_{c}-\alpha_{2}K(x)$ and let $I = (-\infty, \omega)$ be the interval where $g(x)<0$: then any non-zero solution of (8.14) has at most one zero in $I$. Indeed, let $w(x)$ be solution of (8.14) and $x_{0}$ one of its zeros. Put

Taking into account Eq (8.14), we have

So $W$ is non-decreasing on $I$. If $W$ has another zero $x_{1}\in I$, then $W(x) = 0$ on $[x_{0}, x_{1}]$. Thus $w$ is constant on $[x_{0}, x_{1}]$, namely $w = 0$ (from (8.14)) on $[x_{0}, x_{1}]$ and consequently on $I$. The contradiction shows that $w$ has at most one zero in I.

Since $w(x)$ converges to zero when $x\rightarrow\infty$, there exists $x^{*}$ large enough such that we can assume that Eq (8.14) becomes:

Let $c>2\sqrt{D}$. Since $g(x)\rightarrow 1$ at infinity, there exist two linearly independent solutions of Eq (8.17) given by ^[42]

where $\lambda_{1}$ and $\lambda_{2}$ are solution of the algebraic equation

If $\lambda_{1}>\lambda_{2}$, then the general solution of (8.17) can be written as

where $k_{1}$ and $k_{2}$ are real constants. We can choose $k_{1}<0$ and $k_{2}>0$ such that $w(x)$ has an only zero $x_{0}$ and $w(x)<0$ for $x<x_{0}$, $w(x)>0$ for $x>x_{0}$. In addition, $w(x)$ behaves as $v_{c}(x)$ if $c> 2\sqrt{D}$: so that, we can have the estimation $w(x)\leq P(x, 0)\leq v_{c}(x)$ for the initial condition.

If $c = 2\sqrt{D}$, then $\lambda_{1} = \lambda_{2}$. Then the qualitative behavior of the solutions of Eq (8.17) as $x\rightarrow\infty$ is determined by

with $k_{1}$ and $k_{2}$ real constants. We can choose $k_{1}<0$ and $k_{2}>0$ such that $w(x)$ has an only zero $x_{0}$ and $w(x)<0$ for $x<x_{0}$, $w(x)>0$ for $x>x_{0}$. This result can be achieved choosing both the constants $k_{1}$ and $k_{2}$ positive. In addition, $w(x)$ behaves as $v_{c}(x)$ if $c = 2\sqrt{D}$: so that, we can have the estimation $w(x)\leq P(x, 0)\leq v_{c}(x)$ for the initial condition.

Lemma 6. Let $z_{1}(x) = \max(0, w_{c}(x))$, and $z_{2}(x) = v_{c}(x)$. If

then the solution of the Cauchy problem for Eq (8.4) with the initial condition $P^{0}(x)$ satisfies the estimate

for all $t>0$.

Theorem 7. There exist positive GTW solutions of Eq (8.3) for all $c\geq 2\sqrt{D}$. Positive GTW converging to zero as $x\rightarrow\infty$ do not exist for $c<2\sqrt{D}$.

Proof. The existence of GTWs for all $c\geq 2\sqrt{D}$ follows from the previous lemma. Indeed, consider solution of Eq (8.3) in the form $P(x, t) = P(x-ct, t)$. Then

It follows from Lemma 3 that there exists an $\omega$-limit solution $P_{c}(x, t)$ of Eq (8.20) such that

for all $t\in\mathbb{R}$. In order to construct this solution, consider the solution $P(x, t)$ of Eq (8.20) with an initial condition $P^{0}(x)$ which satisfies the inequality $z_{1}(x)\leq P^{0}(x)\leq z_{2}(x)$ for all $x$. Let $t_{n}\rightarrow\infty$ as $n\rightarrow\infty$. Consider next solutions $P_{n}(x, t)$ with the initial conditions $P_{n}^{ 0} = P(x, t_{n})$. Obviously, each of them is defined for $t>-t_{n}$. A locally convergent subsequence of the sequence of functions $P_{n}(x, t)$ is a solution of Eq (8.20) defined for all $t\in\mathbb{R}$. It satisfies inequality (8.21). It can be easily verified that it is a GTW with the speed $c$.

Suppose now that there exists a positive GTW $P_{c}(x, t)$, converging to $0$ as $x\rightarrow\infty$, with a speed $c<2\sqrt{D}$. Then $P_{c}(x-ct, t)$ satisfies Eq (8.20). Let us take $c<c_{0}<2\sqrt{D}$ and consider the equation

It has a solution $P_{0}(x)$ that is non-monotone and unstable ^[41]. Moreover, when $x$ tends to infinity, then $P_{0}(x)\sim exp(-c_{0}x/2)\sin(ax)$, where $a = \sqrt{D|c_{0}^{2}/4-1|}$. Therefore, equation

has a solution $P^{*}(x, t) = \varepsilon P_{0}(x-(c_{0}-c)t)$, where $\varepsilon$ is a positive constant. Let $x = N_{1}$ and $x = N_{2}$ be two consecutive zeros of the function $P_{0}(x)$ such that $P_{0}(x)$ is positive between them. Then $P^{*}(x, t)$ is a solution of the boundary value problem for Eq (8.23) in the domain

with the zero boundary conditions. For $\varepsilon$ small enough, similarly to (8.9) we can obtain the inequality

If $m_{a}(t)$ is the maximal solution of the equation

then

Since $c_{0}>c$ and $P_{c}(x, t)$ converges to zero as $x\rightarrow\infty$, then the last inequality contradicts the assumption that $P_{c}(x, t)$ is a GTW with the speed $c$.

This contradiction proves the theorem.

8.2.1.   Numerical simulations

In Figure 5 we present the results of numerical simulations of Eq (2.10) in one space dimensions. The kernel is the function $\phi_4(x)$. If the support of $\phi$ is sufficiently small, then there is a classical traveling wave propagating with a constant speed. Figure 5(a) shows the solution $P(x, t)$ of Eq (2.10) with the initial condition which has a bounded support. The solution represents two waves propagating in the opposite directions. It is interesting to note that the wave is not monotone with respect to $x$. If we increase the support of the function $\phi$, then the homogeneous in space stationary solution $P_2$ loses its stability and a periodic in space structure appears. In this case we observe propagation of a periodic wave (Figure 5(b)-(c)). Figure 5(d) shows the solution $P(x, t)$ with the exponentially decaying initial conditions. Therefore, depending on the values of parameters equation (2.10) can have solutions of different types. For all other parameters fixed, usual traveling waves (with a constant speed and profile) are observed for sufficiently small values of $h$. Periodic traveling waves exist for sufficiently large $h$. Transition from simple to periodic waves occurs due to the essential spectrum crossing the imaginary axis. The stationary solution homogeneous in space $P_2$ loses its stability resulting in the appearance of a stationary periodic solution. The traveling wave connects the constant value $P_0$ for $x = + \infty$ with this periodic solution for $x\longrightarrow -\infty$ . The waves in (a)-(b)-(c) (in which the initial data has a bounded support) move at their minimal speed $c = 2\sqrt{D} = 2$ and this speed does not depends on $h$; while, the speed of the wave (d) is greater than the minimal speed $c>2\sqrt{D}$ and depends on $h$. Though we consider in numerical simulations a finite interval, if it is sufficiently large, then the solution can approach the corresponding GTW.

9.   The role of the memory of past information

Here we assume that the information is not only spatially and but also temporally non–local: the subjects take their decisions not only on the information on the current state of the system but also on the past states. In other words, we include the memory of the subjects concerning the past information about vaccine side-effects. We consider in this section the 1D unbounded domain, i.e., $\Omega = \mathbb{R}$.

These assumptions yield:

with $W(\tau)$ a delaying kernel, i.e., a positive function such that $\int_{0}^{+\infty}W(\tau)d\tau = 1$. Using scaling (2.9), we obtain:

with $t>0$, $x\in \Omega$.

Here we will focus on the so called acquisition–fading kernel (AFK) ^[46]:

with $0<b<d$. This memory kernel, which is such that $W(0) = 0$ (i.e., absence of information of the current state of the process), models the process of temporal acquisition of information (with a rapid timescale $1/d$) followed by a fading of the memory (with a 'slow' timescale $1/b$).

By applying the linear chain trick (see the Appendix) one gets the following equivalent system:

9.1.   Stability of stationary solutions

Model (9.4) is remindful of systems generating complex self-organized patterns in biochemistry and in population biology ^[22,47], since only one of the involved state densities is endowed of diffusion. The aim of this section is to investigate the possible bifurcation from the homogenous equilibria of (9.4) of such self-organized structures.

It is easy to verify that model (9.4) admits the following homogeneous in space stationary solutions

Also here we will focus on $\mathbf{E}_2$. Considering the case where $\Omega$ is unbounded and linearizing (9.4) at $\mathbf{E}_{2}$ yields the following eigenvalue problem:

whose essential spectrum is given by:

Thus, one obtains the following $\xi$-family of $\lambda$ polynomials

where

Thus $\mathbf{E}_{2}$ is stable iff $a_{0}(\xi)>0$ and

that is:

We recall that a change of sign of $H_1$ in non-spatial systems is associated to the onset of a Hopf bifurcation, a change of sign in a spatiotemporal context give rise to a Hopf-Turing bifurcation ^[48] and finally the only change of sign of $a_0(\xi)$ would give rise to a 'pure' Turing bifurcation ^[48,49].

In this context it is important to distinguish instabilities that genuinely involve the spatio–temporal model form those that would have been also found in the non–spatial model (NSM). As regards the NSM, we observe that

Thus, for Routh-Hurwitz criterion to be satisfied by the NSM, it must be:

Thus:

ⅰ) if $b+d > 1-1/\alpha$, then $H_1^{NSM}>0$ and the equilibrium of the the NSM is LAS;

ⅱ) if $b+d<1-1/\alpha$ then $H_1^{NSM} < 0$ and the equilibrium of the NSM is unstable via Hopf bifurcation.

The above result on $H_1^{NSM}$ actually implies that $H_1(\xi)>0$. Indeed from (9.11), we have that $H_1(\xi)$ is composed by three addenda: $(d+b)D^2 \xi^4 >0$, $(b+d)^2 D\xi^2 >0$ and

So no Turing-Hopf bifurcation can be generated if the equilibrium of the NSM is LAS.

Thus, apparently the only route to instability at $\mathbf{E}_2$ is via changes of sign of $a_0(\xi)$. However it is easy matter to verify that the condition $a_0(\xi)>0$ is nothing else that the LAS condition of the intermediate equilibrium $P_2 = 1/\alpha$ defined in (5.2) of our model (2.10) without memory.

Summarizing, we may state that the introduction of temporal non-locality via the AFK does not make $\mathbf{E}_2$ unstable unless:

ⅰ) $b+d \ge 1-1/\alpha$, i.e., the NSM is unstabilized at $E_2$ by the memory effect;

ⅱ) there exist $\xi$ such that $a_0(\xi)<0$, i.e., the original model without memory is unstable.

The above analysis, with simple changes, can also be repeated in the case of local and of global information.

In the Figure 6 we give an illustration of the stability character of the state of equilibrium $\mathbf{E}_2$ for the system (9.4), with $\phi(x) = \phi_4(x)$, choosing a small perturbation as initial data. In Figure 6(a) we see that the solutions evolve in time to reach the stable equilibrium $\mathbf{E}_2$, while in Figure 6(b) the stability is lost and we see the convergence to a periodic spatial structure. This depends on the diffusion coefficient $D$, the ratio $\rho = \alpha_2/\alpha_1$ and the parameter $h$.

Even if the presence of GTW has not been proved analytically for the system (9.4), in Figure 7 we provide an example of possible GTW in the case of system (9.4). The kernel $\phi(x)$ is the function $\phi_{4}(x)$ and the initial data decays exponentially. We focus our attention on delay parameters $b$ and $d$. If $b$ and $d$ are small enough, then spatial periodic structures appear, while increasing the value of $b$ and $d$, this kind of structure disappear.

9.2.   Traveling wave with local information

Now, we describe wave solutions for Eq (9.4), setting the Dirac delta $\delta(x)$ as kernel, i.e., solutions of the following system

Let $\textbf{E} = (P, Z, M)^T$ be the variable vector of system (9.12) such that

with $\mathbf{E}_{-}$ unstable and $\mathbf{E}_{+}$ stable equilibrium state. After substitution of the wave variable $z = x-ct$ in model (9.12), we obtain the following system of ordinary differential equations:

where $P(x-ct) = p(z)$, $Z(x-ct) = \omega(z)$ and $M(x-ct) = m(z)$. When $\alpha<1$, recall that $\textbf{E}_0$ is unstable and $\textbf{E}_1$ is stable with respect to the system (9.4). Thus, the following result holds ^[44]:

Theorem 8. If $0<\alpha<1$, there exists a traveling wave solution of velocity $c$ $\forall\ c\in[c^{(1)}, +\infty[$, connecting equilibrium $\mathbf{E}_{1}$ and equilibrium $\mathbf{E}_{0}$, i.e., a vector function $\mathbf{E}(x-ct)$, solution of system (9.14) on the real line $]-\infty, +\infty[$ and satisfying Eq (9.13) with $\mathbf{E}_{-} = \mathbf{E}_{1}$ and $\mathbf{E}_{+} = \mathbf{E}_{0}$.

If $\alpha>1$, then equilibrium $\mathbf{E}_1$ becomes unstable while equilibrium $\textbf{E}_2$ is stable with respect to model (9.4). Therefore, it can affirmed that ^[44]:

Theorem 9. If $\alpha>1$, there exists a traveling wave solution of velocity $c$ $\forall\ c\in[c^{(1)}, +\infty[$, connecting equilibrium $\mathbf{E}_{2}$ and equilibrium $\mathbf{E}_{0}$, i.e., a vector function $\mathbf{E}(x-ct)$, solution of system (9.14) on the real line $]-\infty, +\infty[$ and satisfying Eq (9.13) with $\mathbf{E}_{-} = \mathbf{E}_{2}$ and $\mathbf{E}_{+} = \mathbf{E}_{0}$; moreover, traveling wave solution between $\mathbf{E}_{1}$ to $\mathbf{E}_{0}$ cannot exist.

In Figure 8, we show traveling wave solution for system (9.4) with minimal speed $c = c^{(1)} = 2\sqrt{D}$, connecting the equilibrium $\mathbf{E}_2$ at $-\infty$ and equilibrium $\mathbf{E}_0$ at $+\infty$.

10.   Vaccine awareness campaigns: Stability Analysis of Homogeneous Steady States

A relevant task of a Public Health System (PHS) is that of investing to continuously sustain the vaccine uptake in the population. In the scenario we considered, namely prevention of an infection which has been eliminated, this is particularly important. Indeed, in such a case parents are scarcely motivated to vaccinate due to the low perceived risk, which is mirrored by the fact that $\vartheta(M_i) = \vartheta_0$.

The action of the PHS in enacting vaccine awareness campaigns was modelled in ^[14] by amending the basic imitation dynamics of strategies. Briefly, they assumed that the action by the PHS allows a steadily positive flux $\gamma A$ switching per time unit from the Non-vaccinator to the Vaccinator strategy.

Extending the above ideas ^[5,14,15,31] to the present spatiotemporal setting yields the model:

where $\gamma$ is a positive constant.

Keeping expressions (2.3) and (2.4), using $A = 1-P$ and scaling as (2.9), from (10.1) we obtain:

with $t>0$, $x\in\Omega$ and Neumann condition (2.11).

A first effect of public action can immediately be seen: the NVE $P_0 = 0$ is no more an equilibrium point.

Moreover, system (10.2) has two spatially homogenous equilibria, namely the PVE $P_1 = 1$ and the equilibrium $P_{\gamma} \in (1/\alpha, 1)$ where $P_{\gamma}$ is the solution of:

i.e.,

It is easy to verify that:

ⅰ) if $0 < \gamma < (\alpha -1)$ then $P_{\gamma} \in (1/\alpha, 1);$ i.e., $P_{\gamma}$ is epidemiologically meaningful;

ⅱ) it is $d/d\gamma P_{\gamma}>0$;

ⅲ) it is $P_{\gamma} = P_2 = 1/\alpha$, if $\gamma = 0$. The previous results show that the new behavioral equilibrium $P_{\gamma}$, which replaces the 'old' equilibrium $P_{2}$, allows a higher vaccine uptake than $P_{2}$ thanks to the univocal effect of public intervention.

11.   A preliminary simulation of the impact of Public Health System Action

Before investigating the effects of the awareness intervention by the PHS on the stability of the equilibria of (10.2) (and its dependence on the distinct hypotheses on the spatial information kernel), here we illustrate, by means of a simulation, how a Public Health System intervention to favour vaccination can impact on the behavior of the system. We assume that subjects make their decisions by taking into the account their memory. In particular we choose the scenario depicted in the in the panel (a) of Figure 7 a GTW whose features are deeply influenced by the temporal non–locality. We assume that initially there is no public health intervention, which starts at about two third of the simulation interval and grows up to a maximal value $\gamma_{Max}$:

Thus, we simulate the equation:

where $W(\tau)$ is the delay function defined in (9.3) and $\gamma(t)$ the function defined in (11.1). Figures 9 and 10 (view from the top) show that initially the GTW pattern coexists with a rapidly increasing uniform pattern. This pattern relatively soon destroys the GTW and, any case, it immediately increases the minimum value attained by the GTW.

11.1.   Local information

Under purely local information (i.e., $\phi(x) = \delta(x)$), Eq (10.2) becomes:

Proceeding as in the case without vaccine awareness campaign, the following global stability results hold:

Theorem 10. If $\alpha >1$ and $\gamma \in (0, \alpha -1)$ then $P_{\gamma}$ is globally stable (and $P_1$ is unstable).

If $\alpha >1$ and $\gamma \ge \alpha -1$ then $P_1$ is globally stable.

In the trivial case where $0 \le \alpha \le 1$ then $P_1$ is globally stable independently of $\gamma \ge 0$.

The previous results straightforwardly extend to the present spatially structured case the findings in ^[14], showing that the behavioral equilibrium $P_{\gamma}$ is always GAS whenever it exists, while when it disappears the PVE inherits its global stability. As for possible heteroclinic connections between $P_{\gamma}$ and $P_1$, introducing the traveling variable $z = x-ct$, the following negative result is established:

Theorem 11. There are no monotone decreasing traveling wave solutions for Eq (11.3) connecting equilibrium $P_1$ at equilibrium $P_{\gamma}$ with positive speed.

11.2.   Non–local information

11.2.1.   The case of $\Omega$ unbounded

As noted above, assuming that $L\gg 1$, we have the following equation

with $x\in\Omega = \mathbb{R}$, $t>0$ and Neumann condition (2.11).

The associated linearized equation at a generic homogeneous steady state $P^{*}$ reads

and the related eigenvalue problem has the form:

If $P^{*} = P_{\gamma}$, then we obtain the following eigenvalue

To determine the sign of

it is useful to take into the account that

yielding

i.e.,

It is easy to show that if $\gamma<\alpha -1$ then $\chi<0$. Thus, eigenvalues (11.7) become

Therefore, we have demonstrated the following results:

ⅰ) since $\chi_{\gamma}<0$, for both $\phi(x) = \phi_1(x) = C_g e^{-ax^2}$ and $\phi(x) = \phi_2(x) = C_e e^{-a|x|}$ it is

so that $P_{\gamma}$ is LAS;

ⅱ) for $\phi(\xi) = \phi_3(x) = C_N (1-(|x|/N)^2) Hev(N-|x|)$ and $\phi(\xi) = \phi_4(x) = (1/2h)Hev(h-|x|)$ the equilibrium $P_{\gamma}$ can be unstable.

Moreover, if $P^{*} = P_{1} = 1$, consider the eigenvalue problem (11.6), then we obtain

If $\gamma>\alpha -1$, then $\nu(\xi)$ is negative for any $\xi$ and equilibrium $P_{1}$ is LAS.

11.2.2.   The case of bounded $\Omega$

Assuming that the set $\Omega$ is bounded, we obtain the following equation from (10.2):

and Neumann conditions (2.11).

We will consider here $\Omega = [-L, L]$.

The eigenvalues problem associated at the linearization about a generic steady state $P^*$ is given by

Since for the equilibrium $P_{1} = 1$ the nonlocal term disappears, the ensuing local stability analisys is straightforward and it is omitted. Thus we will only focus on the equilibrium $P_{\gamma}$.

Assuming that the kernel $\phi(x)$ satisfy conditions (6.21) and (2.7) i.e., function $\phi$ is periodic with period $2L$, even and normalized, if $P^{*} = P_{\gamma}$, we look for solution of the following eigenvalue problem

in the form

Then the boundary conditions are satisfied. Taking into account that the function $\phi$ is periodic and even, the m–th eigenvalue reads as follows:

where $\phi_{m}$ is defined in (6.23). We note that $\nu_{m}$ is negative for $m = 0$ and for $m$ sufficiently large.

If $\phi_{m}$ is a positive function, then the corresponding eigenvalue $\nu_{m}$ is also positive. If the kernel $\phi$ has a bounded support, then, depending on the parameters, the corresponding eigenvalues can be negative and stability is lost.

11.3.   Global information

In case of global information,

with $x\in\Omega = \mathbb{R}$, $t>0$ and Neumann condition (2.11), it is trivial to show that $P_{\gamma}$ is LAS. Indeed, setting $W(x) = \cos(\xi_m x)$ as eigenfunction of the problem:

we get the following eigenvalues:

12.   Comparison with the Theory of Diffusion of the innovations

In this section we aim at highlighting similarities and key differences between our models (2.1) and (10.1) and the Theory of Innovation Diffusion (TID) by Mahajan and others ^[34,32]. For the sake of the notation simplicity we mainly consider non-spatial models.

The TID models focus on the dynamics of the adoption of innovation. Defining $Y(t)$ as the fraction of adopters of the innovation at time $t$, and $U(t) = 1-Y(t)$ the fraction of 'non–adopters' the basic family of models is the following:

where $g(t)>0$, which was initially interpreted as an 'hazard of adoption' ^[36], has later been compared by Capasso and Zonno to the force of infection of epidemic models ^[33], with whom it has striking analogies. Thus, we will call it 'Force of Innovation Adoption' (FIA).

The positivity of the FIA $g(t)$ has an important consequence for our comparison: all specific models belonging to the family (12.1) are such that $Y(t)\rightarrow 1^{-}$, i.e., there is a unique equilibrium point, which can be termed 'all adopters'. As a consequence, assuming $Y = P$, all possible analogies with our family of models is limited to the unrealistic cases where the 'all vaccinators' equilibrium is GAS. For example, the family of models (12.1) for $g(t) = \gamma$ corresponds to our model in the trivial case where there is no contagion of ideas. The positivity of $g(t)$ and its interpretation as a 'force' of infection also stress the key difference between the family of models we investigated here and the family of model (12.1): in the scenario we investigated there is a double contagion of ideas, i.e., a double flux: ⅰ) from $A$ to $P$; ⅱ) from $P$ to $A$. On the contrary, the family (12.1) implies an uni–directional contagion of ideas. This is a direct consequence of the deeply different nature of the underlying 'social processes'. However, it is of interest to investigate formal analogies between the two theories. The specific instances of imitation game-like equations of TID ^[32,34] are obtained phenomenologically by assuming that $g(t)$ is an analytical function of $Y$ ^[34]:

where $a_k \ge 0$ can also be function of time ^[33,35]. For example, assuming that the non–linear terms are null, $g(t) = a_0+a_1 Y(t)$, yields the most popular popular model of the TID: the Bass model ^[34,36], which reads as follows:

The non–negativity of the coefficients $a_k$ would preclude formal analogies with our models. However, maintaining the constraint of positivity of the FIA $g(t)>0$ but relaxing for $k\ge 2$ the hypothesis that $a_k \ge 0$ yields the following $g(t)$ that allows a formal analogy with some models of the family of models considered in this work:

with, however, the following sharp constraint ($\alpha(Y)$ is non–decreasing):

that again leads to the global stability of the 'all adopters' ('all vaccinators') model.

The important case of linear–affine $\alpha(y)$ can be expressed within the frame of the TDI family of models as follows:

Note that this generalized framework is such that: ⅰ) can be framed with the epidemiological theory by Capasso where forces of infection can also be non–monotone ^[6], which can also be applied to diffusion of information ^[33]; ⅱ) it remains a formal model of a uni-dimensional flow.

Namely, a FIA $g(t)$ of the type (12.4) can be read in terms of TID as an initially beneficial effect of the product qualities, possibly followed, for larger $Y$, by a partial mitigation of the 'enthusiasm' for the product due to the spread of information of its possible (real of presumed) defects. Thus, a reduced FIA is observed.

As far as the space is concerned, let us first consider the Fisher-Kolgomorov–like model (3.1), for which the equilibrium 'all vaccinators' is GAS. If we further assume $\alpha_0 = 0$ then we have a unidirectional flux from the group 'non-vaccinators' to the group 'vaccinators'. In such a case our model (3.1) is equivalent to the Mahajan spatiotemporal model of Innovation Diffusion ^[34].

More interestingly, in ^[35,33] Capasso proposed, in the framework of an SIR-like diffusion of innovation model, a family of non-local models of the FIA. This family of models is non-locally depending on the adopters, again in the context of a uni-directional flux. In our case, instead, we have two fluxes and it is the force of infections from the adopters towards non–adopters that is non–local.

13.   Discussion and concluding remarks

The role played by the spatial structure of information used by vaccination decision makers is a main topic of behavioral epidemiology. In this work we have investigated, within a spatial framework based on classical diffusion, the effects of three different structures of information ((ⅰ) purely 'local' information, 'local' plus a 'global', country-wide, average information, (ⅲ) a mix of local and non-local information) on the dynamics of vaccine uptake in absence of the infection. As a consequence, given the background of low incentive to immunization, the dynamics of VAEs emerge as the key determinants of vaccination decisions and collective coverage. This possibly represents a main case of the current public health landscape in modern industrialized countries, where a number of vaccine preventable infectious disease were successfully eliminated 'locally' (think e.g., to polio which was locally eliminated long ago in most Europe), but infection re-emergence must be prevented by means of large vaccine coverage.

We focused our analyses on the pattern and properties - namely stability, bifurcation, existence of classical and generalized traveling waves, effects of temporal delays, and the effects of awareness campaigns by the Public Health System - of the key space-homogeneous equilibrium solutions, that we termed the 'behavioral equilibria'.

Our main results were as follows.

As regards the stability properties of the behavioral equilibrium, these show a nice interplay between the form of the spatial information kernel and the nature of the spatial domain (bounded vs unbounded). Under “purely local information” the BE is always LAS (and GAS) whenever it exists. Under 'local and global' information the BE is always LAS whenever it exists. However, convergence is slower than in the case of purely local information. Under 'mixed' information results are more articulated. For unbounded domains stability will prevails independently of the shape of the spatial information kernel under large values of the diffusion coefficients. In particular kernels that are strictly decreasing in the distance from the local position are always stabilising independently of the spatial kernel. However, bounded-support kernels (e.g., strictly positive up to some threshold distance and zero thereafter) can yield instability at low levels of the diffusion coefficient when the strength of non-local information tend to prevail on local one.

Interestingly, the instability caused by the presence of non–local information can generate generalized traveling waves characterised by more or less pronounced oscillations as well as other not static spatial patterns. The onset of these traveling waves depends on the interplay between behavioral parameters and the structure of the spatial kernel. This means that the fraction of individuals favourable to vaccination can show oscillations. Then we investigated the case where agents also use past - and not only current - information about VAEs in forming their perceptions of risk. We showed that in the realistic case where this temporal non-locality is modelled by the acquisition–fading kernel no further specific instability can arise w.r.t. the non–spatial scenario.

However, the presence of information memories can remarkably impact on the 'shape' of the generalised traveling waves. Indeed, for certain combination of the characteristic times of the memory kernel, the spatial oscillations internal to the propagating front are so that, at each time, there are large zones where the proportion favorable to vaccination is close to zero. This can compromise the herd immunity of the population.

Remarkably, instabilities and emergence of generalised traveling waves are cleared out soon after the start of vaccine awareness campaigns enacted by the public health system, in line with the intuitions supplied in ^[14].

We briefly mention that the above illustrated mathematical behaviors could also be read, from theoretical physics viewpoint, as phase changes and phase transitions (PTs), an increasingly important concept in ME ^[50,51]. Indeed, on the one hand in general the kinetic of PTs is most often characterized by traveling waves ^[52,53,54]. On the other hand, all imitation games (as our model) where individual switch between two mutually exclusive strategies have some analogy with the well-known Ising model ^[55]. In our case, one can 'read' $P(x, t) \approx 1$ as a phase rich in pro–vaccine subjects (an 'ordered' phase) and $P(x, t) \approx 0$ as a phase poor in pro–vaccine subjects (a 'disordered phase').

Finally, we highlighted fundamental differences and some analogies between the phenomena in study and those investigated by the Theory of Innovation Diffusion. Basically in classical TID models, there is a uni–directional contagion of ideas between 'non-adopters' and 'adopters'. As a consequence all classical models of TID predicts global stability of the full adoption of the innovation. This correspond to the pathological case where $P = 1$ is Globally Asymptotically Stable.

This work obviously has a number of limitations.

A key limitation of this study is the hypothesis that at the macroscale the agents' mobility can be approximated by classical diffusion i.e., by random walk around geographic space. This is a coarse approximation, thus our model is essentially a benchmark for providing clear-cut baseline results. They will have to be validated against more robust hypotheses. Indeed, real patterns of mobility of human individuals are complex ^[56] and can introduce remarkable implications for infection patterns, such as unexpected phase transitions ^[57]. A promising recent route is to model human mobility as a superdiffusive process ^{[58,59,60,61]}. Superdiffusion allows a potentially more realistic and flexible representation of human spatial mobility while, at the same time, remaining analytically tractable. Thus superdiffusion is a natural candidate for future extensions of the present work. In particular, Brockmann and Hufnagel ^[62] have modelled a double chemical reaction leading to equations similar to a classical imitation game in case of superdiffusive mobility of molecules and we plan to follow their approach in a follow-up work.

A second limitation of this study lies in the specific case-study, namely a vaccination scenario in the absence of the infection. As mentioned in the introduction, this case is a central. However, it is equally important the scenario where the disease is continuously re–introduced. This, for example, is the case of measles, which has caused sizable epidemics in recent years ^[8,9,63]. Further forthcoming work will therefore be devoted to the study of the spatio-temporal interplay between the vaccine opinion dynamics and infection spread.

Acknowledgments

We thank the four anonymous referees that with their important comments helped us to significantly improve this manuscript. The publication has been prepared with the support of the RUDN University Program 5-100, and the French-Russian project PRC2307. The publication was also supported by LYSM-INDAM and GNFM-INdAM.

Conflict of interest

The authors declare there is no conflict of interest.

APPENDIX: reduction to two differential equations of the information model with memory

Here, it is shown that the integro-differential information model (9.2) is equivalent to the differential system (9.4) with the delaying kernel $W(\tau)$ defined in (9.3). First, consider the function

which is clearly the solution of the following linear differential equation:

Given $W(\tau)$ defined in (9.3), consider

i.e.,

Therefore,

i.e.,