# A geometric analysis of the impact of large but finite switching rates on vaccination evolutionary games

Rossella Della Marca, Alberto D'onofrio, Mattia Sensi, Sara Sottile

<figure><img src="image_1.png" /></figure>

► To cite this version:

Rossella Della Marca, Alberto D'onofrio, Mattia Sensi, Sara Sottile. A geometric analysis of the impact of large but finite switching rates on vaccination evolutionary games. Nonlinear Analysis: Real World Applications, 2023, 75, pp.103986. 10.1016/j.nonrwa.2023.103986. hal-04086364

HAL Id: hal-04086364

https://inria.hal.science/hal-04086364

Submitted on 2 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A geometric analysis of the impact of large but finite
switching rates on vaccination evolutionary games

Rossella Della Marca¹, Alberto d’Onofrio²,*, Mattia Sensi³,⁴,*, Sara Sottile⁵

¹Mathematics Area, SISSA – International School for Advanced Studies,
via Bonomea 265, I-34136 Trieste, Italy
rossella.dellamarca@sissa.it

²Department of Mathematics and Geosciences, University of Trieste,
via Alfonso Valerio 12/1, 34127 Trieste, Italy
alberto.donofrio@units.it (*corresponding author)

³MathNeuro Team, Inria at Université Côte d'Azur, 2004 Rte des Lucioles, 06410 Biot, France,

⁴Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino Italy
mattia.sensi@polito.it (*corresponding author)

⁵Department of Mathematics, University of Trento, Via Sommarive 14, 38123 Povo - Trento, Italy
sara.sottile@unitn.it

Abstract

In contemporary society, social networks accelerate decision dynamics causing a rapid switch of opinions in a number of fields, including the prevention of infectious diseases by means of vaccines. This means that opinion dynamics can nowadays be much faster than the spread of epidemics. Hence, we propose a Susceptible-Infectious-Removed epidemic model coupled with an evolutionary vaccination game embedding the public health system efforts to increase vaccine uptake. This results in a global system “epidemic model + evolutionary game”. The epidemiological novelty of this work is that we assume that the switching to the strategy “pro vaccine” depends on the incidence of the disease. As a consequence of the above-mentioned accelerated decisions, the dynamics of the system acts on two different scales: a fast scale for the vaccine decisions and a slower scale for the spread of the disease. Another, and more methodological, element of novelty is that we apply Geometrical Singular Perturbation Theory (GSPT) to such a two-scale model and we then compare the geometric analysis with the Quasi-Steady-State Approximation (QSSA) approach, showing a criticality in the latter. Later, we apply the GSPT approach to the disease prevalence-based model already studied in (Della Marca and d’Onofrio, *Comm Nonl Sci Num Sim*, 2021) via the QSSA approach by considering medium-large values of the strategy switching parameter.
# 1 Introduction

The increasing spread of hesitancy and refusal of vaccines is a major challenge for global public health. This problem was originally born in the field of prevention of childhood diseases and of influenza but it has largely been observed also during the current pandemic of COVID-19 [42]. Focusing on childhood diseases vaccines, we may say that one of major determinants of this problem is the phenomenon of “pseudo-rational” exemption to vaccination. This phenomenon consists in the fact that parents both overweight real or supposed side effects of vaccines and underweight real disease-related risks [33, 43]. This causes the public health systems (PHSs) to spend considerable energies and budget to mitigate the impact of this phenomenon, typically by means of vaccine awareness public campaigns aimed at increasing the vaccine uptake.

In consequence of the outbreak of exemption to vaccination, a new scientific discipline has been developed: the Behavioural Epidemiology of Infectious Diseases (BEID) [33,43], whose aim is the inclusion in epidemic models of the description of human decision making (concerning, e.g., vaccination choices, social distancing, mobility patterns).

Not surprisingly, Game Theory has an important place in the context of BEID. Namely, the dynamics of vaccine decision is very frequently modelled as an imitation evolutionary game [4, 17] which can be represented as an infection of ideas process [43]. Of course, also the modeling of the above-mentioned efforts of PHSs aimed at increasing the vaccine uptake have been introduced in the BEID literature [18].

In the papers [4, 17, 18] the implicit focus was on slow changes of vaccine strategy. In the age of social media this implicit assumption is often unrealistic. Indeed, social networks accelerate decision dynamics causing a rapid switch of opinions in a number of fields, including the prevention of infectious diseases, for example in vaccination campaign. This means that opinion dynamics can nowadays be faster than the spread of epidemics. We are witnessing what has been defined as an “exponential growth in public opinion channels” [29], which leads public opinion to be extremely volatile on many key subjects [3], such as politics [3, 27] and vaccines [10, 36]. This scenario, unimaginable until few years ago, must guide contemporary models in BEID.

Recently, Della Marca and d'Onofrio [14] explored the impact of the above-mentioned volatility on the modelling of public response to vaccine awareness campaigns for favouring vaccine uptake. Since the evolutionary vaccination game is endowed of a parameter that tunes the velocity of strategy change, under the hypothesis that such parameter is extremely large, they applied a Quasi-Steady-State Approximation (QSSA) of the model [18]. This resulted in a Susceptible–Infectious–Removed (SIR) epidemic model with a nonlinear dependence on the control: the PHS effort to increase the vaccine uptake. This control was designed via optimal control approach and numerically implemented also via heuristic global optimization methods.

In the present work, we drastically depart from the paper [14] in a number of modeling and methodological points. Firstly, here we hypothesize that the switching from the vaccine refusal to vaccine acceptance is influenced by the information on current disease incidence. This is a more realistic hypothesis since the most widely diffused information on the spread of an infectious disease
is not the prevalence but the incidence. Secondly, apart from the extreme cases that the velocity of strategy change is small or is practically infinite, we consider an intermediate case: the study of the impact of a large but finite switching velocity via the Geometric Singular Perturbation Theory (GSPT) [22].

GSPT [22] is a powerful approach to model phenomena evolving on multiple time scales [5, 23, 26, 31]. Such separation in time scales differing by many orders of magnitude is quite common in real world scenarios (e.g., chemical oscillations, neuroscience [15, 16, 37, 40, 41], electromechanical devices [13], lasers, ecology, celestial mechanics, pattern formation [31]). In particular, it has been applied to epidemics models in which immunity windows and demographic turnover are much longer than infectious periods [2, 12, 24, 25], and to discrete time epidemic models [6–9].

In this work, in line with what was done e.g. in the paper [39], we take into the account that the time scale of vaccination strategy changes, although very fast, nonetheless is not instantaneous. Hence, we model the volatility of strategy switching by means of rigorous GSPT. We compare our geometric analysis to a classical QSSA, highlighting when the two approaches lead to the same conclusion, and when the QSSA instead fails to reproduce the exact dynamics for a large but finite strategy switching rate. This failure, as we will illustrate, is due to a delayed loss of stability of the critical manifold in our model. We exploit the so-called entry-exit function [32, 34, 35, 38] to characterize this crucial part of the dynamics of our model. Moreover, in the final part of this work we briefly apply the GSPT also to the model heuristically inferred in the paper [18].

The manuscript is organized as follows. In Section 2, we introduce an evolutionary vaccination game where the information on the spread of the disease concerns the disease incidence; in Section 3, we perform a qualitative analysis of the model in the case of low/medium strategy switching rate; in Section 4, we consider the case of very large (infinite) switching rate and apply the QSSA; in Section 5, we investigate the case of large but finite switching rate by using the GSPT approach and compare the results to those obtained by the QSSA; in Section 6, we provide some numerical simulations; in Section 7, we apply the GSPT approach to the evolutionary vaccination game studied in the paper [18], where the information on the spread of the disease concerns the disease prevalence; we conclude in Section 8, summarizing our main results and providing inspiration for future research.

## 2 Background of Evolutionary Vaccination Game as a process of Mutual Infection of ideas

Let us consider the SIR-like model describing the dynamics of a vaccine-preventable endemic childhood disease under voluntary vaccination choices [17, 18]:

$$ 
\begin{aligned}
\dot{S} &= \mu(1 - p(\tau)) - \mu S - \beta SI, \\
\dot{I} &= \beta SI - (\nu + \mu)I,
\end{aligned}
 $$

where: i) $\tau$ denotes the (slow) time variable; ii) $S$ and $I$ represent the fraction of susceptible and infectious individuals within the population at time $\tau$; iii) birth and natural mortality rates are equal to a value $\mu$; iv) $p$ is the time-dependent vaccine uptake of newborns; v) $\beta$ is the disease transmission rate; vi) $\nu$ is the recovery rate.
Note that (iii) implies a stationary population, so that we can neglect the equation ruling the dynamics of the fraction of removed individuals: $R = 1 - S - I$. Observe that, independently of the dynamics of $p$, the basic reproduction number is

$$ \mathcal{R}_0 = \frac{\beta}{\nu + \mu}. \qquad (1) $$

To model the impact of human decision making on the vaccination choices, we assume that the population of parents is proportional to the total (constant) population and is divided into two groups [4, 17, 18]: “pro-vaccine” and “anti-vaccine”. The first group is given by parents who are in favour of vaccines and vaccinate their children ($p$); the second one is given by parents who are hesitant or overtly against vaccination and, as a consequence, do not vaccinate their children ($a = 1 - p$).

The evolution of $p$ follows an imitation game dynamics that could be inferred by employing an economics-oriented approach based on the concept of payoffs [4, 17] or, in alternative, a statistical physics-oriented approach [11, 43]. We follow the second approach that we consider to be much clearer than the first one in the present context. The basic concept is that the dynamics of $p$ and $a$ are ruled by a “double contagion” of ideas between the two involved groups. This approach yields the following family of models [43]:

$$ \begin{aligned} \dot{p} &= k_1 \bar{\theta}(M_d) pa - k_1 \bar{\alpha}(M_v) ap, \\ \dot{a} &= -k_1 \bar{\theta}(M_d) pa + k_1 \bar{\alpha}(M_v) ap, \end{aligned} $$

where: i) the “force of infection” concerning the switch from the strategy “anti-vaccine” to the strategy “pro-vaccine” is $k_1 \bar{\theta}(M_d)p$, where $M_d(\tau)$ is an information variable on the extent of the disease status in the community (e.g., incidence, prevalence); ii) the “force of infection” concerning the switch from the strategy “pro-vaccine” to the strategy “anti-vaccine” is $k_1 \bar{\alpha}(M_v)a$, where $M_v(\tau)$ is an information/“rumors” variable on the extent of vaccine-related side effects; iii) the parameter $k_1$ is a time scale tuning parameter that characterizes the velocity of strategy switching in the population of parents. In the following, we assume that both $\bar{\theta}(\cdot)$ and $\bar{\alpha}(\cdot)$ are linear-affine functions, namely we set:

$$ \bar{\theta}(M_d) = \theta_0 + \theta_1 M_d, \quad \bar{\alpha}(M_v) = \alpha_0 + \alpha_1 M_v, \qquad (2) $$

with $\theta_0 \ge \alpha_0$. The action of the PHS to favour the vaccine uptake is simply modeled as an additional switch, say $k_1 \gamma_1 a$, from the strategy “anti-vaccine” to the strategy “pro-vaccine”. This provides the model:

$$ \begin{aligned} \dot{p} &= k_1(\theta_0 + \theta_1 M_d)pa - k_1(\alpha_0 + \alpha_1 M_v)ap + k_1\gamma_1 a, \\ \dot{a} &= -k_1(\theta_0 + \theta_1 M_d)pa + k_1(\alpha_0 + \alpha_1 M_v)ap - k_1\gamma_1 a. \end{aligned} $$

Taking into the account that $a = 1 - p$ yields the following imitation game equation for $p$:

$$ \dot{p} = k_1 p(1-p)(\theta_0 - \alpha_0 + \theta_1 M_d - \alpha_1 M_v) + k_1 \gamma_1 (1-p). \qquad (3) $$

It is easy to see that one cannot identify all of these parameters:

$(k_1, \alpha_1, \gamma_1).$

The reason is simple: $k_1$ never appears “alone” so, unless one has an *a priori* knowledge of the value of $\gamma_1$ or of $\alpha_1$, $k_1$ cannot be obtained from a parameter estimation of epidemiological data using our model.
**Remark 1.** The scenario is similar to the Malthusian parameters b and m:

$$ \dot{x} = (b - m)x = rx. $$

Unless one has an external measure of b or of m, one cannot identify both b and m from data concerning $x(t)$. One can only fit the difference

$$ r = b - m. $$

Thus, the best option is to rewrite the imitation game equation (3) rescaling the parameters by $\alpha_1$, as follows

$$ k = k_1\alpha_1, \quad \delta = \frac{\theta_0 - \alpha_0}{\alpha_1}, \quad \theta = \frac{\theta_1}{\alpha_1}, \quad \gamma = \frac{\gamma_1}{\alpha_1}. \qquad (4) $$

As mentioned in Section 1, we focus here on the relevant case in which the main information on the disease available to the population is the number of new cases. Thus, $M_d$ is a measure of the available information on the disease incidence at time $\tau$. In the case in which the decisions are taken by only considering current information, we obtain

$$ M_d(\tau) = \beta S(\tau) I(\tau). \qquad (5) $$

We further assume that the information on vaccine side-effects is proportional to the vaccine uptake of newborns, namely $M_v(\tau) = p(\tau)$, like in the papers [14, 17, 18]. We thus obtain the following complete model:

$$ \dot{S} = \mu(1 - p) - \mu S - \beta SI, \qquad (6a) $$

$$ \dot{I} = \beta SI - (\nu + \mu)I, \qquad (6b) $$

$$ \dot{p} = k(p(1-p)(\delta + \theta\beta SI - p) + \gamma(1-p)). \qquad (6c) $$

## 3 An Evolutionary Game in which strategy switching depends on the incidence of the disease

Let us start by investigating the case of low/medium rate of strategy switching ($k \gg 1$), corresponding to classical societies with far less volatile opinions.

The model (6) presents two disease-free equilibria. The first one is a disease-free state in which all the newborns are vaccinated:

$$ E^P = (0, 0, 1). $$

It is shown in Appendix A.1 that high values of $\gamma$, i.e. $\gamma > \bar{\gamma}$, where

$$ \bar{\gamma} = 1 - \delta, $$

ensure the global attractivity of $E^P$. Conversely, when $\gamma < \bar{\gamma}$, $E^P$ is unstable.
The second disease-free equilibrium is

$$E^0 = (1 - p^0, 0, p^0),$$

where

$$p^0 = \frac{\delta + \sqrt{\delta^2 + 4\gamma}}{2}, \qquad (7)$$

which exists only when $p^0 < 1$, namely $\gamma < \bar{\gamma}$. This disease-free state bifurcates from the equilibrium $E^0$ at $\gamma = \bar{\gamma}$.

Let us introduce the threshold value

$$\gamma_c = p_c(p_c - \delta) < \bar{\gamma}, \qquad (8)$$

with

$$p_c = 1 - \frac{1}{\mathcal{R}_0},$$

and $\mathcal{R}_0$ as given in (1). It is possible to show that if $\gamma > \gamma_c$ (i.e. $p^0 > p_c$), then $E^0$ is globally attractive (see Appendix A.2). Further, if $\gamma < \gamma_c$ (i.e. $p^0 < p_c$), then $E^0$ is unstable. When $E^0$ becomes unstable, a unique (and epidemiologically meaningful) endemic equilibrium

$$E^* = (S^*, I^*, p^*)$$

appears by a transcritical bifurcation at $\gamma = \gamma_c$. The components of $E^*$ read

$$S^* = \frac{1}{\mathcal{R}_0}, \quad I^* = \frac{\mu}{\nu + \mu} (p_c - p^*), \quad p^* = \frac{\delta + \theta \mu p_c + \sqrt{(\delta + \theta \mu p_c)^2 + 4\gamma(1 + \theta \mu)}}{2(1 + \theta \mu)}. \qquad (9)$$

Straightforward calculations show that $I^* > 0$ is equivalent to $p^0 < p_c$, that is $\gamma < \gamma_c$. The local stability of the equilibrium $E^*$ may vary with the model parameters, in particular Hopf bifurcations can occur, as we will prove shortly in Theorem 2. First, let us introduce, for simplicity, the notation

$$q_1 = \mu + \beta I^*, \quad q_2 = (\nu + \mu)\beta I^*, \quad q_3 = \frac{\gamma}{p^*} + p^*, \qquad (10)$$

with $I^*$ and $p^*$ as given in (9).

**Theorem 2.** The endemic equilibrium $E^*$ of model (6) is locally asymptotically stable (LAS) if $W \ge 0$, with

$$W = q_1^2 q_3 + \theta \mu p^* \beta I^* (\beta I^* - \nu) + 2 \sqrt{q_1 q_2 q_3 (q_1 q_3 + \theta \mu p^* \beta I^*)},$$

and $q_1, q_2, q_3$ as given in (10). Otherwise, if $W < 0$, then there exist two values $k_1, k_2$ with $0 < k_1 < k_2$, such that $E^*$ is unstable for $k \in (k_1, k_2)$, whereas it is LAS for $k < k_1$ or $k > k_2$. Hopf bifurcations occur at $k = k_i, i = 1, 2$. Moreover, in such a case, if $k \in (k_1, k_2)$ then the orbits are oscillatory in the sense of Yakubovich.
Proof. The Jacobian matrix of system (6) evaluated at the equilibrium $E^*$ reads

$$J = \begin{pmatrix} -q_1 & -(\nu + \mu) & -\mu \\ \beta I^* & 0 & 0 \\ k p^*(1-p^*)\theta\beta I^* & k p^*(1-p^*)\theta(\nu + \mu) & -k(1-p^*)q_3 \end{pmatrix},$$

leading to the characteristic polynomial

$$p(\lambda) = \lambda^3 + a_1\lambda^2 + a_2\lambda + a_3,$$

with

$$a_1 = k(1-p^*)q_3 + q_1 > 0, \quad a_2 = k(1-p^*)(q_1q_3 + \theta\mu p^*\beta I^*) + q_2 > 0, \quad a_3 = k(1-p^*)q_2(q_3 + \theta\mu p^*) > 0,$$

where $q_1, q_2, q_3$ are given in (10). The positiveness of the coefficients of $p(\lambda)$ rules out, by Descartes rule of signs, the possibility of real positive eigenvalues, so that stability losses of the endemic state can only occur via Hopf bifurcations. More precisely, according to Routh-Hurwitz theorem, $E^*$ is LAS if and only if $a_1a_2 - a_3$ is positive, equivalently written as

$$f(k) = b_0k^2 + b_1k + b_2 > 0,$$

with

$$b_0 = (1 - p^*)^2 q_3 (q_1 q_3 + \theta \mu p^* \beta I^*) > 0, \quad b_1 = (1 - p^*)(q_1^2 q_3 + \theta \mu p^* \beta I^* (\beta I^* - \nu)), \quad b_2 = q_1 q_2 > 0.$$

We chose the strategy switching rate, $k$, as bifurcation parameter because it affects the stability but not the existence of $E^*$.

Thus, if $b_1 \ge 0$ and/or if $\Delta = b_1^2 - 4b_0b_2 \le 0$, then $f(k) > 0$ independently of $k$.

Otherwise, if $b_1 < 0$ and $\Delta > 0$, then $f(k)$ has two positive roots $k_1, k_2$, with $0 < k_1 < k_2$, such that $f(k) < 0$ for $k \in (k_1, k_2)$ and $f(k) > 0$ for $k < k_1$ or $k > k_2$. In such a case, at $k_i$, $i = 1, 2$, the test for non-zero speed is fulfilled:

$$f'(k)|_{k=k_i} = \pm\sqrt{\Delta} \neq 0.$$

Note that, by simple algebra, one can write

$$\Delta = (b_1 - 2\sqrt{b_0b_2})(b_1 + 2\sqrt{b_0b_2}).$$

It follows that, if

$$b_1 + 2\sqrt{b_0b_2} < 0,$$

then $b_1 < 0$ and $\Delta > 0$. Vice versa, if

$$b_1 + 2\sqrt{b_0b_2} \ge 0,$$

then $b_1 \ge 0$ and/or $\Delta \le 0$.

Finally, as regards the Yakubovich oscillatority, it holds that: i) the orbits of the system are bounded; ii) the endemic equilibrium $E^*$ is unstable when $f(k) < 0$; iii) the disease-free equilibria are unstable. Thus, we may apply the Yakubovich theorem [20, 21]. The claim follows. $\square$
<figure><img src="image_3.png" /><figcaption>Figure 1: Hopf bifurcation locus for the equilibrium $E^*$ of model (6). Roots $k_1, k_2$ of the Routh-Hurwitz function $f(k)$ in Theorem 2 as functions of the factor of perceived risk of infection $\theta \in [0, 25000/(\nu + \mu)]$. Parameter values as given in Table 1, with $\psi = 0.3$.</figcaption></figure>

A bifurcation diagram showing the Hopf bifurcation locus $k_1, k_2$ as functions of the factor of perceived risk of infection $\theta$ is given in Fig. 1.

Regarding the Yakubovich oscillatority, this intuitively means that for sufficiently large time all the state variables are permanently oscillating, with regular or irregular oscillations (periodic, quasi-periodic or chaotic). Formally, it holds that [20, 21]

$$-\infty < \liminf_{\tau \to +\infty} X(\tau) < \limsup_{\tau \to +\infty} X(\tau) < +\infty,$$

for $X \in \{S, I, p\}$. Note that this is a global result, unlike the Hopf bifurcation theorem, which is local.

# 4 The case of very large k: Quasi-Steady-State Approximation

Let us investigate now the case of very large rate of strategy switching, corresponding to the extreme case that the velocity of switching is practically infinite.

We consider system (6) with

$$k = \frac{1}{\varepsilon},$$
and $\varepsilon > 0$ very small. The ensuing system is a slow-fast system with two slow variables, $S$ and $I$, and one fast variable, $p$.

If we assume that $\varepsilon \ll 1$, then a QSSA for $p$ can be used, which yields

$$p(1-p)(\delta + \theta\beta SI - p) + \gamma(1-p) = \varepsilon\dot{p} \approx 0.$$

As a consequence, in the limit $\varepsilon \to 0$, $p$ is the solution of the following algebraic equation:

$$0 = (1-p)(\gamma + (\delta + \theta\beta SI)p - p^2), \quad (11)$$

to be solved under the constraint $0 \le p \le 1$. Equation (11) has two solutions: $p = 1$, and the unique positive solution of

$$\gamma + (\delta + \theta\beta SI)p - p^2 = 0.$$

Summarizing, $p$ tends to

$$p = \zeta(S, I) = \min \left( 1, \frac{\delta + \theta\beta SI + \sqrt{(\delta + \theta\beta SI)^2 + 4\gamma}}{2} \right) \quad (12)$$

because it is stable, contrary to $p = 1$ which is unstable. In other words:

i) if $p = \zeta$, then $\dot{p} = 0 \quad \forall \varepsilon > 0$;

ii) if $p \in [0, \zeta)$, then

$$\lim_{\varepsilon \to 0} \dot{p} \to +\infty;$$

iii) if $p \in (\zeta, 1)$, then

$$\lim_{\varepsilon \to 0} \dot{p} \to -\infty.$$

Thus, model (6) reduces to the following bidimensional model:

$$\begin{aligned} \dot{S} &= \mu(1 - \zeta(S, I)) - \mu S - \beta SI, \\ \dot{I} &= \beta SI - (\nu + \mu)I, \end{aligned} \quad (13)$$

with $\zeta$ defined in (12).

Model (13) admits an unique disease-free equilibrium

$$E^0 = (1 - p^0, 0),$$

with $p^0 = \zeta(S, 0)$. When $\gamma < \gamma_c$, with $\gamma_c$ as given in (8), it also admits the endemic equilibrium

$$E^* = \left( \frac{1}{\mathcal{R}_0}, I^* \right),$$

where $I^*$ is the unique positive solution of

$$\zeta\left(\frac{1}{\mathcal{R}_0}, I^*\right) = p_c - \frac{\nu + \mu}{\mu} I^*. \quad (14)$$
Indeed, the l.h.s. of (14) is an increasing function of $I^*$, the r.h.s. of (14) is a decreasing function of $I^*$, and $\zeta(1/\mathbb{R}_0, 0) = p^0 < p_c$. Note that, with a slight abuse of notation, we denote by $E^0$ and $E^*$ also the corresponding equilibria of the model (6) in the case $\varepsilon \ll 1$ (see Section 3).

By using arguments similar to those of paper [19], it can be shown that, if $\gamma > \gamma_c$ then $E^0$ is globally asymptotically stable (GAS), see Appendix A.3; instead, if $\gamma < \gamma_c$ then $E^0$ is unstable and $E^*$ is GAS in the positively invariant region

$$ \Omega^* = \{(S, I) | S \ge 0, I > 0, S + I \le 1, S \le 1 - p^0\} $$

(see Appendix A.4).

The proposed QSSA provides a model of the spread and control of an SIR-like infectious disease that extends the one in the paper [19] to the important case where the information on the disease spread is the incidence, not the prevalence. As in the paper [19], oscillations do not occur since the rate of strategy change is too large. The physical reason of this lack of limit cycles and other oscillating structures is that here we are in a regime of extremely volatile public opinion. In other words, there is no opinion-induced delay with respect to the information on the disease spread.

**Remark 3.** The above procedure can easily be adapted to a far more general case in which $\bar{\theta}(\cdot)$ and $\bar{\alpha}(\cdot)$ are nonlinear, in place of (2). This would result in a more general nonlinear relationship between $p(t)$, the instantaneous incidence ($\beta SI$) and the control parameter $\gamma$: $p = F(\beta SI; \gamma) \in [0, 1]$.

# 5 An approach based on Geometric Singular Perturbation Theory

In this section, we focus on a different approach to the model under analysis, through the use of techniques from GSPT.

We denote by $\tau$ the slow time variable and with $t = k\tau = \tau/\varepsilon$ the fast time variable. In a more classical GSPT notation, system (6) can be rewritten, in the fast time scale $t$, as

$$ \begin{aligned} S' &= \varepsilon(\mu(1-p) - \mu S - \beta SI), \\ I' &= \varepsilon(\beta SI - (\nu + \mu)I), \\ p' &= p(1-p)(\delta + \theta\beta SI - p) + \gamma(1-p). \end{aligned} \quad (15) $$

Note that we denote by $\dot{X}$ the derivative of the variable $X$ with respect to $\tau$ and by $X'$ the derivative with respect to $t$, where $X \in \{S, I, p\}$.

**Remark 4.** The unit of measure of $k$ is $1/\text{time}$. This implies that, when we apply a change in the time coordinate, bringing system (6) to system (15), the resulting $t = k\tau = \tau/\varepsilon$ is dimensionless. Numerical simulations of (15) should be handled carefully, since $t$ does not have a time dimension; thus, we need to rescale the time accordingly.

The results obtained with GSPT are asymptotic as $\varepsilon \to 0$. Specifically, this means that for each result there exists a $\varepsilon_0 > 0$, which is often not explicitly quantifiable, such that the result holds for $0 < \varepsilon < \varepsilon_0$. This translates, considering $k = 1/\varepsilon$, to a validity of our results for $k \in (1/\varepsilon_0, +\infty)$.
In some cases, the QSSA is in perfect agreement with the results obtained with GSPT. In other cases, namely for orbits which pass exponentially close in $\epsilon$ to the manifold $\{p = 1\}$, the QSSA is not able to replicate the delayed loss of stability of that same manifold. This behaviour is clearly visible in Fig. 2, in which different situations have been explored.

The critical manifold is given by

$$ \begin{aligned} C_0 &= \{S \ge 0, I \ge 0, S + I \le 1, p \in [0, 1] \mid p(1-p)(\delta + \theta\beta SI - p) + \gamma(1-p) = 0\} \\ &= \{p = 1\} \cup \{p^2 - (\delta + \theta\beta SI)p - \gamma = 0\} =: A \cup B. \end{aligned} \quad (16) $$

We first make the expression of B explicit. To remain in the biologically feasible region, we impose $p \in [0, 1]$, from which it follows that $p = \zeta(S, I)$, as given in (12). Also, it must be

$$ SI \le \frac{1 - \delta - \gamma}{\theta\beta}, $$

where we have assumed that $2 - \delta - \theta\beta SI \ge 0$ on an open interval in $(0, 1)^2$. Moreover, we need that $1 - \delta > \gamma$, mirroring the necessary conditions coming from the paper [18].

Computing the Jacobian of the fast system (15) with $\epsilon = 0$, we notice that there are two zero eigenvalues, corresponding to the slow variables S and I. The third eigenvalue is

$$ \lambda = (p^2 - (\delta + \theta\beta SI)p - \gamma) + (1-p)(\delta + \theta\beta SI - 2p). $$

On B, the eigenvalue is

$$ \lambda = (1-p)(\delta + \theta\beta SI - 2p). $$

Note that $1 - p \ge 0$ and, on B, $\delta + \theta\beta SI - 2p = -\sqrt{(\delta + \theta\beta SI)^2 + 4\gamma} < 0$, thus $\lambda < 0$ and B is always locally attractive when it exists, independently of the value of $\mathcal{R}_0$.

Let us now focus on the behaviour of the system when $p \approx 1$ (i.e. close to A). If we sum the first two equations in (6), we obtain

$$ \dot{S} + \dot{I} \le -\mu(S + I), $$

thus $S + I$ converges to zero. We remark that on A the susceptible population S can only decrease, since this set represents the situation in which there are no vaccine sceptical people, whereas I is not always decreasing, although the sum $S + I$ is. The eigenvalue on A is

$$ \lambda = p^2 - (\delta + \theta\beta SI)p - \gamma. $$

Thus,

$$ \left\{ \begin{array}{ll} \lambda > 0 & \text{if } \frac{\delta + \theta\beta SI + \sqrt{(\delta + \theta\beta SI)^2 + 4\gamma}}{2} < p \le 1, \\ \lambda < 0 & \text{if } 0 \le p < \frac{\delta + \theta\beta SI + \sqrt{(\delta + \theta\beta SI)^2 + 4\gamma}}{2}. \end{array} \right. $$

Since the corresponding eigenvalue $\lambda$ changes its sign, the region A is attractive until the intersection with the curve (12), it then becomes repelling and, after a delay, the dynamics lands on the curve B. Depending on the value of $\mathcal{R}_0$ and $S_{in}$, where $S_{in}$ indicates the entrance of an orbit in a neighbourhood of $\{p = 1\}$, we can observe two different behaviours:
* Case I: if $\mathcal{R}_0 < 1$, or $\mathcal{R}_0 > 1$ and $S_{in} < 1/\mathcal{R}_0$, then $I(\tau)$ is decreasing and the orbits converge to $B$ after a delay $T_E$ determined with the entry-exit function $\int_0^{T_E} \lambda(\tau)d\tau = 0$; hence the exit time $T_E$ is given implicitly by

$$ (1 - \delta - \gamma)T_E = \theta\beta \int_0^{T_E} S(\tau)I(\tau)d\tau. $$

A similar entry-exit phenomenon, with orbits eventually landing on a different branch of the critical manifold, was already observed in the paper [1]. We note that, however, we have no explicit formula for $S(\tau)I(\tau)$ on $p = 1$. There is, to the best of the authors' knowledge, no theoretical result that justifies the canard-like behaviour of orbits remaining in a neighbourhood of an unstable branch in this setting. However, from the lower dimensional case of Section 7, for which the delayed loss of stability is theoretically foreseeable, and from our simulations, we conjecture that this is always the case, for the parameter values we are interested in;

* Case II: if $\mathcal{R}_0 > 1$ and $S_{in} > 1/\mathcal{R}_0$, then $I(\tau)$ changes its monotonicity and the convergence towards $B$ only happens after a short excursion away from it. In this case, the known formulas for entry-exit functions cannot be applied.

# 6 Numerical simulations: the impact of finite $k$

For numerical simulations, we use epidemiological parameters compatible with a vaccine-preventable endemic childhood disease [14, 17, 18]. In particular, we keep fixed the parameters values: $\nu = 52$ years$^{-1}$, $\mu = 1/78$ years$^{-1}$, $\mathcal{R}_0 = 18$, $\varepsilon = 0.1/\nu$. At variance, we vary the values of the imitation game parameters $\delta$, $\gamma$ and $\theta$, that are chosen in the following way:

* since the value of $\delta$ affects the value of $\gamma_c$, recall (8), we define it as $\delta = \psi p_c$, where $\psi \in \{0.3, 0.8\}$;

* we then compute the value of $\gamma_c$ and we set $\gamma = \gamma_c/2$;

* the most complex case concerns finding the value of $\theta$. We cannot use the same value used in the paper [14], say $\theta_{pre} = 450$, since in that case the perceived risk of infection is proportional to $I(\tau)$, whereas here it is proportional to $\beta S(\tau)I(t)$. We impose that, if $S(\tau)$ is close to its endemic equilibrium value $S^* = 1/\mathcal{R}_0$ and the value of $I(\tau)$ is equal to that in the prevalence-based case, then the two “risks of infection” must be equal each other, namely $\theta\beta S^*I(\tau) = \theta_{pre}I(\tau)$. This gives in turn a baseline value for $\theta$ in the present case, say $\theta_{bas} := \theta_{pre}/(\mu + \nu)$. We eventually set $\theta = \phi\theta_{bas}$, where $\phi \in \{10, 50\}$.

As regards the initial conditions of the state variables, we assume that they are at the endemic equilibrium of the model (6) with constant vaccination $p(\tau) \equiv 0.75$; we choose this value because it is large enough, and it is representative of an initial condition in which three quarters of parents used to vaccinate their children. Therefore, $S(0) = 1/\mathcal{R}_0$, $I(0) = \mu(p_c - p(0))/(\mu + \nu)$, $p(0) = 0.75$.

The above-mentioned values of the parameters and the initial conditions are summarized in Table 1.
<table><caption>Table 1: Parameters values and initial conditions used for numerical simulations of models (13) and (15).</caption><thead><tr><td>Parameter</td><td>Formula</td><td>Value</td></tr></thead><tbody><tr><td>ν</td><td>-</td><td>52 years<sup>-1</sup></td></tr><tr><td>μ</td><td>-</td><td>1/78 years<sup>-1</sup></td></tr><tr><td>ε</td><td>0.1/ν</td><td>0.0019 years</td></tr><tr><td>k</td><td>1/ε</td><td>520 years<sup>-1</sup></td></tr><tr><td>β</td><td>R<sub>0</sub>(μ + ν)</td><td>936.23 years<sup>-1</sup></td></tr><tr><td>p<sub>c</sub></td><td>1 - 1/R<sub>0</sub></td><td>0.944</td></tr><tr><td>R<sub>0</sub></td><td>-</td><td>18</td></tr><tr><td>ψ</td><td>-</td><td>{0.3; 0.8}</td></tr><tr><td>δ</td><td>ψp<sub>c</sub></td><td>{0.283; 0.756}</td></tr><tr><td>θ<sub>pre</sub></td><td>-</td><td>450</td></tr><tr><td>θ<sub>bas</sub></td><td>θ<sub>pre</sub>/(μ + ν)</td><td>8.652 years</td></tr><tr><td>φ</td><td>-</td><td>{10; 50}</td></tr><tr><td>θ</td><td>φθ<sub>bas</sub></td><td>{86.52; 432.59} years</td></tr><tr><td>γ<sub>c</sub></td><td>p<sub>c</sub>(p<sub>c</sub> - δ)</td><td>{0.624; 0.178}</td></tr><tr><td>γ</td><td>γ<sub>c</sub>/2</td><td>{0.089; 0.312}</td></tr><tr><td>S(0)</td><td>1/R<sub>0</sub></td><td>0.056</td></tr><tr><td>I(0)</td><td>μ(p<sub>c</sub> - p(0))/(μ + ν)</td><td>4.793 · 10<sup>-5</sup></td></tr><tr><td>p(0)</td><td>-</td><td>0.75</td></tr></tbody></table>

Since $I(\tau)$ becomes exponentially small in $\varepsilon$, meaning $I = e^{-K/\varepsilon}$ for some $K > 0$, we apply a change of variables to system (15), as reported in Appendix A.5, in order to reduce its numerical stiffness. Moreover, recall from Remark 4 that the time should be rescaled.

In Figs. 2-4 we display the numerical solutions of models (13) and (15), as well as the critical manifold (16), for the four possible combinations of the factors $\psi$ and $\phi$ defining the imitation game parameters $\delta$ and $\theta$, as indicated in Table 1. Namely: (a) $\psi = 0.3$, $\phi = 10$; (b) $\psi = 0.8$, $\phi = 10$; (c) $\psi = 0.3$, $\phi = 50$; (d) $\psi = 0.8$, $\phi = 50$. Specifically, we combine cases in which: i) $\delta$ is considerably smaller than $p_c$ ($\psi = 0.3$) or quite close to it ($\psi = 0.8$); ii) $\theta$ is medium ($\phi = 10$) or relatively large ($\phi = 50$). Note that, in the case of absence of action enacted by the PHS ($\gamma = 0$), $\theta$ is the slope of the reactivity with respect to the information on the disease spread.

We report in Fig. 2 the dynamics of the vaccine uptake of newborns, $p$; in Fig. 3 the dynamics of the fraction of infectious individuals normalized with respect to the endemic value in (9), $I/I^*$; in Fig. 4 the dynamics of the disease incidence, $\beta SI$.

From Figs. 2-4, we note a generalized and very remarkable discrepancy between the solutions by the QSSA (model (13), red lines) and those by the GSPT (model (15), black lines). In any case, the solutions of both the models converge towards the endemic equilibrium. However, in the case (a), the damped oscillations of $p$ predicted by the QSSA have much larger period than those predicted by the GSPT (Fig. 2a). In the cases (b)-(c)-(d), the orbits of the model (15) exhibit transitory oscillations by approaching $\{p = 1\}$, while those of model (13) approach the endemic equilibrium after a short passage near $\{p = 1\}$.
Vaccine uptake of newborns $p$

<figure><img src="image_4.png" /><figcaption>Figure 2: Impact of the finiteness of the rate of switching strategy $k$. Vaccine uptake of newborns as predicted by the QSSA model (13) (red lines) and by the GSPT model (15) (black lines). The black dashed lines represent the critical manifold (16). Panel (a): $\psi = 0.3, \phi = 10$. Panel (b): $\psi = 0.8, \phi = 10$. Panel (c): $\psi = 0.3, \phi = 50$. Panel (d): $\psi = 0.8, \phi = 50$. Other parameter values and initial conditions as given in Table 1.</figcaption></figure>

Further, in the case (a), the disease prevalence and incidence predicted by the model (15) resemble the shape of their QSSA counterpart (Fig. 3a and Fig. 4a). At variance, in the cases (b)-(c)-(d), the disease prevalence and incidence predicted by the GSPT showcase effect of the delayed loss of stability of the system through a slow passage near $I = 0$ and $SI = 0$, respectively. From Fig. 3 we also note that for both the models varying $\psi$ and $\phi$ (and hence, $\delta$ and $\theta$) affects the timing of the oscillations of the disease prevalence. In particular, increasing the values of $\psi$ and/or $\phi$ increases the oscillation period.

A less pronounced but nonetheless important discrepancy is also observed in Fig. 2 between the solutions of GSPT (model (15)) and the values of the critical manifold (16) at the solutions of GSPT (black dashed lines). The discrepancy is evident in all the cases except for the case (a), where the orbits of both models (13) and (15) never approach $\{p = 1\}$, and the critical manifold (16) gives an excellent approximation of the behaviour of the model (15). In the cases (b)-(c)-(d), the discrepancy is not limited to very short periods but it extends from some months (Fig. 2b) to many years (Fig. 2d), which can have a remarkable impact from the Public Health viewpoint.
Disease prevalence $I/I^*$

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 3: Dynamics of the fraction of infectious individuals as predicted by the QSSA model (13) (red lines) and by the GSPT model (15) (black lines), and normalized with respect to the endemic equilibrium value $I^*$, given in (9). Panel (a): $\psi = 0.3, \phi = 10$. Panel (b): $\psi = 0.8, \phi = 10$. Panel (c): $\psi = 0.3, \phi = 50$. Panel (d): $\psi = 0.8, \phi = 50$. Other parameter values and initial conditions as given in Table 1.</figcaption></figure>
Disease incidence $\beta SI$

<figure><img src="image_9.png" /><figcaption>Figure 4: Dynamics of the disease incidence $\beta SI$ as predicted by the QSSA model (13) (red lines) and by the GSPT model (15) (black lines). Panel (a): $\psi = 0.3, \phi = 10$. Panel (b): $\psi = 0.8, \phi = 10$. Panel (c): $\psi = 0.3, \phi = 50$. Panel (d): $\psi = 0.8, \phi = 50$. Other parameter values and initial conditions as given in Table 1.</figcaption></figure>
# 7 Reconsidering the case in which prevalence is the main information

In this section, we reconsider our initial imitation game equation (3) and briefly examine the case that the main information on the spread of the disease is the prevalence at time $\tau$, namely

$$M_d(\tau) = I(\tau).$$

This case was partially investigated in the paper [14] by using a QSSA approach. We now complete the analysis by using tools from GSPT.

Note that we continue to assume that $M_v$ is given by the vaccine uptake of newborns: $M_v(\tau) = p(\tau)$.

The corresponding model in the fast time scale $t$ now reads:

$$ 
\begin{aligned}
S' &= \varepsilon(\mu(1-p) - \mu S - \beta SI), \\
I' &= \varepsilon(\beta SI - (\nu + \mu)I), \\
p' &= p(1-p)(\delta + \theta_{\text{pre}}I - p) + \gamma(1-p),
\end{aligned} 
\quad (17) $$

where

$$ \theta_{\text{pre}} = \frac{\theta_1}{\alpha_1} $$

is the corresponding factor of perceived risk of infection, and $k$, $\delta$ and $\gamma$ as defined in (4).

We recall the expression of the endemic equilibrium $E^* = (S^*, I^*, p^*)$ of model (17), whose components are given by

$$ S^* = \frac{1}{\mathcal{R}_0}, \quad I^* = \frac{\mu}{\nu + \mu} (p_c - p^*), \quad p^* = \frac{\delta + \frac{\theta_{\mu}}{\mu + \nu} p_c + \sqrt{\left(\delta + \frac{\theta_{\mu}}{\mu + \nu} p_c\right)^2 + 4\gamma \left(1 + \frac{\theta_{\mu}}{\mu + \nu}\right)}}{2\left(1 + \frac{\theta_{\mu}}{\mu + \nu}\right)}. \quad (18) $$

The critical manifold is

$$ 
\begin{aligned}
C_0 &= \{S \ge 0, I \ge 0, S + I \le 1, p \in [0, 1] \mid p(1-p)(\delta + \theta_{\text{pre}}I - p) + \gamma(1-p) = 0\} \\
&= \{p = 1\} \cup \{p^2 - (\delta + \theta_{\text{pre}}I)p - \gamma = 0\} =: A \cup B,
\end{aligned} 
\quad (19) $$

where $B$ is described by

$$ p = \xi(I) = \frac{\delta + \theta_{\text{pre}}I + \sqrt{(\delta + \theta_{\text{pre}}I)^2 + 4\gamma}}{2}. \quad (20) $$

With the same argument employed in Section 5, we obtain that the Jacobian of the fast system (17) with $\varepsilon = 0$ has two zero eigenvalues and a third eigenvalue given by

$$ \lambda = (p^2 - (\delta + \theta_{\text{pre}}I)p - \gamma) + (1-p)(\delta + \theta_{\text{pre}}I - 2p). $$
On B, the eigenvalue is

$$ \lambda = (1 - p)(\delta + \theta_{\text{pre}}I - 2p) < 0, $$

thus B is always locally attractive when it exists, independently of the value of $\mathcal{R}_0$.

We now study the behaviour of the system when $p \approx 1$ (i.e. close to A). As in Section 5, we have that the quantity $S + I$ converges to zero.

The eigenvalue on A is

$$ \lambda = p^2 - (\delta + \theta_{\text{pre}}I)p - \gamma. $$

Thus,

$$ \left\{ \begin{array}{ll} \lambda > 0 & \text{if } \frac{\delta + \theta_{\text{pre}}I + \sqrt{(\delta + \theta_{\text{pre}}I)^2 + 4\gamma}}{2} < p \le 1, \\ \lambda < 0 & \text{if } 0 \le p < \frac{\delta + \theta_{\text{pre}}I + \sqrt{(\delta + \theta_{\text{pre}}I)^2 + 4\gamma}}{2}. \end{array} \right. $$

Since the corresponding eigenvalue $\lambda$ changes its sign, the region A is attractive until the intersection with the curve (20), it then becomes repelling and, after a delay, the dynamics lands on the curve B. Depending on the value of $\mathcal{R}_0$ and $S_{\text{in}}$, where $S_{\text{in}}$ indicates the entrance of an orbit in a neighbourhood of $\{p = 1\}$, we can observe two different behaviours:

* Case I: if $\mathcal{R}_0 < 1$, or $\mathcal{R}_0 > 1$ and $S_{\text{in}} < 1/\mathcal{R}_0$, then $I(\tau)$ is decreasing and the orbits converge to B after a delay $T_E$ determined with the entry-exit function $\int_0^{T_E} \lambda(\tau)d\tau = 0$; hence the exit time $T_E$ is given implicitly by

$$ (1 - \delta - \gamma)T_E = \theta_{\text{pre}} \int_0^{T_E} I(\tau)d\tau. $$

A similar entry-exit phenomenon, with orbits eventually landing on a different branch of the critical manifold, was already observed in the paper [1]. The canard-like behaviour is confirmed applying the results in the paper [30] (see also [28]). We note that, however, we have no explicit formula for $I(\tau)$ on $p = 1$;

* Case II: if $\mathcal{R}_0 > 1$ and $S_{\text{in}} > 1/\mathcal{R}_0$, then $I(\tau)$ changes its monotonicity and the convergence towards B only happens after a short excursion away from it. In this case, the known formulas for entry-exit functions cannot be applied.

In Figs. 5-6, we report the numerical solutions of model (17), as well as the critical manifold (19), and compare them with those by the QSSA approach. We display in Fig. 5 the dynamics of the vaccine uptake of newborns, $p$; in Fig. 6 the dynamics of the fraction of infectious individuals normalized with respect to the endemic value in (18), $I/I^*$. The parameters used are the same as in Table 1, with the exception of the value of $\theta_{\text{pre}}$. Indeed, in this case we set

$$ \theta_{\text{pre}} = \phi450, \quad \text{with } \phi = \{10, 50\}. \qquad (21) $$

From Fig. 5, we note that even in the scenario of prevalence-based information, the “exact” dynamics of $p$ (black lines) largely differs from the one predicted by the QSSA method (red lines). Moreover, the solutions of the GSPT model (17) may differ from the corresponding values along the critical
Vaccine uptake of newborns $p$

<figure><img src="image_10.png" /><figcaption>Figure 5: Information based on the disease prevalence: impact of the finiteness of the rate of switching strategy $k$. Vaccine uptake of newborns as predicted by the QSSA model proposed in the paper [14] (red lines) and by the GSPT model (17) (black lines). The black dashed lines represent the critical manifold (19). Panel (a): $\psi = 0.3, \phi = 10$ in (21). Panel (b): $\psi = 0.8, \phi = 10$ in (21). Panel (c): $\psi = 0.3, \phi = 50$ in (21). Panel (d): $\psi = 0.8, \phi = 50$ in (21). Other parameter values and initial conditions as given in Table 1.</figcaption></figure>

manifold (19) (black dashed lines) for time intervals of some years (Fig. 5b and Fig. 5d), which may have a remarkable impact from the Public Health viewpoint.

Similarly to what observed in Section 6, from Fig. 6 we also note that for both the models increasing $\delta$ and/or $\theta_{\text{pre}}$ increases the oscillation period of the disease prevalence.
Disease prevalence $I/I^*$

<figure><img src="image_11.png" /><figcaption>Figure 6: Information based on the disease prevalence: dynamics of the fraction of infectious individuals as predicted by the QSSA model proposed in the paper [14] (red lines) and by the GSPT model (17) (black lines), and normalized with respect to the endemic equilibrium value $I^*$ given in (18). Panel (a): $\psi = 0.3, \phi = 10$ in (21). Panel (b): $\psi = 0.8, \phi = 10$ in (21). Panel (c): $\psi = 0.3, \phi = 50$ in (21). Panel (d): $\psi = 0.8, \phi = 50$ in (21). Other parameter values and initial conditions as given in Table 1.</figcaption></figure>

# 8 Concluding remarks

The aims of this work were twofold. To start, we first considered and studied the impact of the use of disease incidence instead of disease prevalence in the context of the vaccination imitation game. This assumption allows for a more realistic model because the incidence of a disease is a more widely diffuse datum than the prevalence. The study of the proposed model was done both for small and for very large parameter $k$, in the second case by using the QSSA.

It is of interest to note that the QSSA provides a model of the spread and control of an SIR-like infectious disease that extends – on a mechanistic ground – the purely phenomenological model [19] to the important case where the information on the disease spread is the incidence, not the prevalence. As in the paper [19], also here oscillations do not occur. However, here we have a rationale: the lack of oscillations is related to the fact that the rate of strategy change is too large. The physical reason of this lack of limit cycles and other oscillating structures is that here we are in a regime of
extremely volatile public opinion. As a consequence, there is no opinion-induced delay with respect to the information on the disease spread. Another important point to stress is that the QSSA allows to precisely identify the role and the impact of the parameter $\gamma$ (which summarizes the efforts of PHS to convince hesitant people to change their vaccine-related decisions) on the vaccine uptake $p$, which turns out to be a nonlinear function of both the disease incidence and of $\gamma$.

Then, we investigate the impact of the boundedness of the large parameter $k$. This allowed us to compare the analytical approaches of QSSA and GSPT.

The main result obtained by means of the GSPT approach lies in the significant extension of the QSSA approach, which is only valid in the limit as $\varepsilon \to 0$ (equivalently, $k \to +\infty$) that may be excessively sharp. As we remarked in Section 5, instead, GSPT provides results for $k$ large but finite (equivalently $\varepsilon$ small but strictly positive), which is more realistic from a behavioural switching point of view. This is crucial for the system under study, due to the presence of a transversal intersection of two parts of the critical manifolds. Orbits travelling close to this intersection do not always follow the stable branches of the critical manifold; rather, they spend a non-negligible amount of time in the vicinity of a repelling branch. This part of the dynamics can be characterized through the use of the so-called entry-exit function, as we showcased in Sections 5 and 7.

In the case analysed in Section 7, the canard-like behaviour can be explained through the application of known analytical results. The model we focus on for most of the work, i.e. the one in Section 5, shares many similarities to the one in Section 7. However, to the best of the authors' knowledge, there is not general analytical result for the behaviour of multiple time scales system near a transcritical point in such a setting. A higher dimensional result in the spirit of the paper [30] would be precious both from a purely theoretical point of view, and from its application potential.

Of course, this study suffers a number of limitations. We may mainly mention three: i) from the epidemiology viewpoint, the adopted model is a non-spatial and deterministic mean-field model; ii) from the behavioural viewpoint, the model is equipped with a relatively simplistic evolutionary model of opinion change; iii) we have no precise, universal way to quantify the range of validity for $k$ (equivalently, $\varepsilon$). Finally, we explicitly stress that our model is far more apt to describe the control of the spread of a childhood infectious disease than to deal with COVID-19, which would require a far more detailed model of both spread and control.

**Acknowledgements.** The present work was performed under the auspices of the National Groups for Mathematical Physics (GNFM) (members: R.D.M. and A.d'O.) and for Mathematical Analysis, Probability and their Applications (GNAMPA) (member: S.S.) of the Italian National Institute for Advanced Mathematics (INdAM).

M.S. and S.S. were supported by the Italian Ministry for University and Research (MUR) through the PRIN 2020 project "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics" (No. 2020JLWP23).
# A Appendix

## A.1 Global stability of $E^P$ for model (6)

Let us assume that $\gamma > \bar{\gamma}$, where $\bar{\gamma} = 1 - \delta$. From the differential equation (6c) it immediately follows that

$$ \dot{p} \geq k(p(\delta - 1) + \gamma)(1 - p), $$

implying that

$$ \liminf_{t \to +\infty} p \geq \min \left( 1, \frac{\gamma}{1 - \delta} \right). $$

Since $\gamma > \bar{\gamma}$, then $\liminf_{t \to +\infty} p = 1$, thereby showing that $E^P$ is GAS.

Conversely, if $\gamma < \bar{\gamma}$, then the instability of $E^P$ easily follows by linearising the differential equation (6c) at $E^P$.

## A.2 Global stability of $E^0$ for model (6)

Assume now that $\gamma_c < \gamma < \bar{\gamma}$, where $\gamma_c = p_c(p_c - \delta)$. Consider the following inequality

$$ \dot{p} \geq k p(1-p) \left( \delta - p + \frac{\gamma}{p} \right), $$

it follows that

$$ \liminf_{t \to +\infty} p = p^0. $$

Moreover, the above minimum limit implies that for large times

$$ \dot{S} \leq \mu(1 - p^0 - S) $$

and in turn that

$$ \dot{I} \leq \beta I(p_c - p^0). $$

Simple algebra shows that $\gamma > \gamma_c$ is equivalent to $p^0 > p_c$, and the global stability of $E^0$ immediately follows.

Conversely, if $\gamma < \gamma_c$ (i.e. $p^0 < p_c$), then the instability of $E^0$ follows from the linearized equation for the infectious fraction: $\dot{i} = \beta i(p_c - p^0)$.

## A.3 Global stability of $E^0$ for model (13)

By defining $\sigma = S + I$, from system (13) we obtain

$$ \dot{\sigma} = \mu(1 - \zeta(S, I)) - \mu\sigma - \nu I \leq \mu(1 - p^0) - \mu\sigma, $$

it follows that

$$ \limsup_{t \to +\infty} \sigma = 1 - p_0. $$

The above maximum limit implies that for large times

$$ \dot{I} \leq \beta I(p_c - p^0). $$

Hence, if $p^0 > p_c$, then $E^0$ is globally attractive.
## A.4 Global stability of $E^*$ for model (13)

The Jacobian matrix of system (13) evaluated at the equilibrium $E^*$ reads

$$J = \begin{pmatrix} -\mu \frac{\partial \zeta}{\partial S} \left( \frac{1}{\mathcal{R}_0}, I^* \right) - \beta I^* - \mu & -\mu \frac{\partial \zeta}{\partial I} \left( \frac{1}{\mathcal{R}_0}, I^* \right) - (\nu + \mu) \\ \beta I^* & -(\nu + \mu) \end{pmatrix}.$$

Since $\text{tr}J < 0$ and $\text{det}J > 0$, the eigenvalues of $J$ have negative real parts. Further, in $\Omega^*$ there are no closed orbits since

$$\text{div} \left( \frac{1}{I}(\dot{S}, \dot{I}) \right) = -\frac{\mu}{I} \frac{\partial \zeta}{\partial S} - \beta - \frac{\mu}{I} < 0.$$

Thus, by the Poincaré-Bendixon thricotomy it follows that $E^*$ is GAS in $\Omega^*$.

## A.5 Change of variables

In order to reduce the numerical stiffness of the model (15), we employ a change of the state variables of the model. Specifically, we introduce:

$$x := \ln(S), \quad y := \ln(I), \quad z := \ln\left(\frac{p}{1-p}\right),$$

i.e.

$$S = e^x, \quad I = e^y, \quad p = \frac{e^z}{1+e^z}.$$

By deriving the new variables with respect to the time $t$, we obtain

$$\begin{aligned} x' &= \varepsilon \left( \mu \frac{e^{-x}}{1+e^z} - \beta e^y - \mu \right), \\ y' &= \varepsilon (\beta e^x - \nu - \mu), \\ z' &= \delta + \theta \beta e^{x+y} - \frac{e^z}{1+e^z} + \gamma \frac{1+e^z}{e^z}. \end{aligned}$$

The obtained system is remarkably less numerically stiff than the original one.

## References

[1] M. A. Achterberg and M. Sensi. A minimal model for adaptive SIS epidemics. *arXiv preprint* 2209.01224, 2022.

[2] M. Aguiar, B. W. Kooi, A. Pugliese, M. Sensi, and N. Stollenwerk. Time scale separation in the vector borne disease model SIRUV via center manifold analysis. *medRxiv preprint* 2021.04.06.21254992, 2021.

[3] A. I. Alberici and P. Milesi. The influence of the Internet on the psychosocial predictors of collective action. *Journal of Community and Applied Social Psychology*, 23(5):373–388, 2013.
[4] C. T. Bauch. Imitation dynamics predict vaccinating behaviour. *Proceedings of the Royal Society B: Biological Sciences*, 272(1573):1669-1675, 2005.

[5] R. Bertram and J. E. Rubin. Multi-timescale systems and fast-slow analysis. *Mathematical Biosciences*, 287:105-121, 2017.

[6] R. Bravo de la Parra, M. Marvá, E. Sánchez, and L. Sanz. Discrete models of disease and competition. *Discrete Dynamics in Nature and Society*, 5310837, 2017.

[7] R. Bravo de la Parra, M. Marvá, E. Sánchez, and L. Sanz. A discrete predator-prey ecoepidemic model. *Mathematical Modelling of Natural Phenomena*, 12(2):116-132, 2017.

[8] R. Bravo de la Parra and L. Sanz. A discrete model of competing species sharing a parasite. *Discrete and Continuous Dynamical Systems - B*, 25(6):2121, 2020.

[9] R. Bravo de la Parra and L. Sanz. Discrete epidemic models with two time scales. *Advances in Difference Equations*, 2021(1):1-24, 2021.

[10] D. A. Broniatowski, A. M. Jamison, S. Qi, L. AlKulaib, T. Chen, A. Benton, S. C. Quinn, and M. Dredze. Weaponized health communication: Twitter bots and russian trolls amplify the vaccine debate. *American Journal of Public Health*, 108(10):1378-1384, 2018.

[11] B. Buonomo, G. Carbone, and A. d'Onofrio. Effect of seasonality on the dynamics of an imitation-based vaccination model with public health intervention. *Mathematical Biosciences and Engineering*, 15(1):299-321, 2018.

[12] C. Castillo-Chavez, D. Bichara, and B. R. Morin. Perspectives on the role of mobility, behavior, and time scales in the spread of diseases. *Proceedings of the National Academy of Sciences*, 113(51):14582-14588, 2016.

[13] A. Cunha Jr, M. Pereira, R. Avanço, A. M. Tusset, and J. M. Balthazar. On the reduction of nonlinear electromechanical systems. *Meccanica*, 57(10):2679-2697, 2022.

[14] R. Della Marca and A. d'Onofrio. Volatile opinions and optimal control of vaccine awareness campaigns: chaotic behaviour of the forward-backward sweep algorithm vs. heuristic direct optimization. *Communications in Nonlinear Science and Numerical Simulation*, 98:105768, 2021.

[15] M. Desroches and V. Kirk. Spike-adding in a canonical three-time-scale model: superslow explosion and folded-saddle canards. *SIAM Journal on Applied Dynamical Systems*, 17(3):1989-2017, 2018.

[16] M. Desroches, B. Krauskopf, and H. M. Osinga. Mixed-mode oscillations and slow manifolds in the self-coupled FitzHugh-Nagumo system. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 18(1):015107, 2008.
[17] A. d’Onofrio, P. Manfredi, and P. Poletti. The impact of vaccine side effects on the natural history of immunization programmes: an imitation-game approach. *Journal of Theoretical Biology*, 273(1):63–71, 2011.

[18] A. d’Onofrio, P. Manfredi, and P. Poletti. The interplay of public intervention and private choices in determining the outcome of vaccination programmes. *PLoS ONE*, 7(10):e45653, 2012.

[19] A. d’Onofrio, P. Manfredi, and E. Salinelli. Vaccinating behaviour, information, and the dynamics of SIR vaccine preventable diseases. *Theoretical Population Biology*, 71(3):301–317, 2007.

[20] D. V. Efimov and A. L. Fradkov. Yakubovich’s oscillatority of circadian oscillations models. *Mathematical Biosciences*, 216(2):187–191, 2008.

[21] D. V. Efimov and A. L. Fradkov. Oscillatority of nonlinear systems with static feedback. *SIAM Journal on Control and Optimization*, 48(2):618–640, 2009.

[22] N. Fenichel. Geometric singular perturbation theory for ordinary differential equations. *Journal of Differential Equations*, 31(1):53–98, 1979.

[23] G. Hek. Geometric singular perturbation theory in biological practice. *Journal of Mathematical Biology*, 60(3):347–386, 2010.

[24] H. Jardón-Kojakhmetov, C. Kuehn, A. Pugliese, and M. Sensi. A geometric analysis of the SIR, SIRS and SIRWS epidemiological models. *Nonlinear Analysis: Real World Applications*, 58:103220, 2021.

[25] H. Jardón-Kojakhmetov, C. Kuehn, A. Pugliese, and M. Sensi. A geometric analysis of the SIRS epidemiological model on a homogeneous network. *Journal of Mathematical Biology*, 83(4):1–38, 2021.

[26] C. K. R. T. Jones. Geometric singular perturbation theory. In L. Arnold, C. K. R. T. Jones, K. Mishaikow, and G. Raugel, editors, *Dynamical systems, Lectures given at the 2nd Session of the Centro Internazionale Matematico Estivo (C.I.M.E.) held in Montecatini Terme, Italy, June 13 - 22, 1994*, volume 1609 of *Lecture Notes in Mathematics*, pages 44–118. Springer, 1995.

[27] J. T. Jost, P. Barberá, R. Bonneau, M. Langer, M. Metzger, J. Nagler, J. Sterling, and J. A. Tucker. How social media facilitates political protest: Information, motivation, and social networks. *Advances in Political Psychology*, 39(S1):85–118, 2018.

[28] P. Kaklamanos, C. Kuehn, N. Popović, and M. Sensi. Entry-exit functions in fast-slow systems with intersecting eigenvalues. *arXiv preprint 2208.11559*, 2022.

[29] D. H. Kamens. Chapter 9 the construction of public opinion. In *Beyond the Nation-State*, pages 257–279. Emerald Group Publishing Limited, 2015.

[30] M. Krupa and P. Szmolyan. Extending slow manifolds near transcritical and pitchfork singularities. *Nonlinearity*, 14(6):1473, 2001.
[31] C. Kuehn. *Multiple Time Scale Dynamics*. Springer, Cham, 2015.

[32] W. Liu. Exchange lemmas for singular perturbation problems with certain turning points. *Journal of Differential Equations*, 167(1):134–180, 2000.

[33] P. Manfredi and A. d'Onofrio, editors. *Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases*. Springer, New York, 2013.

[34] A. I. Neishtadt. Persistence of stability loss for dynamical bifurcations I. *Differential Equations*, 23:1385–1391, 1987.

[35] A. I. Neishtadt. Persistence of stability loss for dynamical bifurcations II. *Differential Equations*, 24:171–176, 1988.

[36] J. On, H.-A. Park, and T.-M. Song. Sentiment analysis of social media on childhood vaccination: Development of an ontology. *Journal of Medical Internet Research*, 21(6):e13456, Jun 2019.

[37] S. Rodrigues, M. Desroches, M. Krupa, J. M. Cortes, T. J. Sejnowski, and A. B. Ali. Time-coded neurotransmitter release at excitatory and inhibitory synapses. *Proceedings of the National Academy of Sciences*, 113(8):E1108–E1115, 2016.

[38] S. Schecter. Exchange lemmas 2: General exchange lemma. *Journal of Differential Equations*, 245(2):411–441, 2008.

[39] S. Schecter. Geometric singular perturbation theory analysis of an epidemic model with spontaneous human behavioral change. *Journal of Mathematical Biology*, 82(6):1–26, 2021.

[40] M. Sensi, M. Desroches, and S. Rodrigues. Slow-fast dynamics in a neurotransmitter release model: delayed response to a time-dependent input signal. *arXiv preprint arXiv:2302.12852*, 2023.

[41] H. Taher, D. Avitabile, and M. Desroches. Bursting in a next generation neural mass model with synaptic dynamics: a slow-fast approach. *Nonlinear Dynamics*, 108(4):4261–4285, 2022.

[42] G. Troiano and A. Nardi. Vaccine hesitancy in the era of COVID-19. *Public Health*, 194:245–251, 2021.

[43] Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d'Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salathé, and D. Zhao. Statistical physics of vaccination. *Physics Reports*, 664:1–113, 2016.