Oscillations in epidemic models with spread of awareness

Abstract

We study ODE models of epidemic spreading with a preventive behavioral response that is triggered by awareness of the infection. Previous studies of such models have mostly focused on the impact of the response on the initial growth of an outbreak and the existence and location of endemic equilibria. Here we study the question whether this type of response is sufficient to prevent future flare-ups from low endemic levels if awareness is assumed to decay over time. In the ODE context, such flare-ups would translate into sustained oscillations with significant amplitudes. Our results show that such oscillations are ruled out in Susceptible-Aware-Infectious-Susceptible models with a single compartment of aware hosts, but can occur if we consider two distinct compartments of aware hosts who differ in their willingness to alert other susceptible hosts.

1 Introduction

The impact of behavioral responses on the progress of an infectious disease has received growing attention in epidemic modeling in recent years (Bagnoli et al. 2007; Epstein et al. 2008; Ferguson 2007; Funk et al. 2010b; Manfredi and d’Onofrio 2013). Responses arising from risk perception of the disease are, for instance, avoidance behavior implying breaking off connections with infectious acquaintances (social distancing), and preventive behaviors, like handwashing or wearing face masks (Kuo et al. 2011; Lau et al. 2010). In epidemic network models, which take into account the contact network in a population, social avoidance has been modeled by means of several mechanisms of dynamic rewiring (Gross et al. 2006; Juher et al. 2013; Llensa et al. 2014; Risau-Gusmán and Zanette 2009; Schwarzkopf et al. 2010). Other approaches cast preventive actions into classic epidemic models by dividing each of the susceptible/infectious/removed classes between responsive and non-responsive individuals (Funk et al. 2010a; Kiss et al. 2010), or by explicitly considering a new class of individuals who are both susceptible and aware, with a lowered susceptibility to infection relative to unaware hosts. This approach leads to the class of SAIS models (Juher et al. 2015; Sahneh et al. 2012, 2014; Sahneh and Scoglio 2011).

Following Funk et al. (2009, 2010a), among others, here we will consider awareness as a state of knowledge about the prevalence of the disease that will both induce a behavioral response in the given host and that this host is willing to transmit to other hosts. To distinguish this more stringent notion of awareness from the one made in the SAIS models of the previous paragraph, we call the latter type of models reactive SAIS models. These models incorporate direct transmission of information about the disease, so that the spread of awareness is somewhat similar but not entirely analogous to the spread of an infectious disease. See Sect. 2 for details.

We are interested in two aspects of the question under what circumstances a behavioral response that is induced by awareness can be an effective control measure: whether such models would predict an elevated epidemic threshold and whether the response would prevent future flare-ups from low endemic levels. Some recent empirical findings appear to show such flare-ups when awareness or at the least adoption of the corresponding behavioral response diminishes. Several public health agencies have reported an increase in some sexual transmitted infections (STIs), particularly among MSM (men who have sex with men) in high-income countries, since the mid-1990s. For instance, the Public Health Agency of Canada reports that syphilis infections are increasing in Canada and, between 2003 and 2012, rates increased by 100% (Totten et al. 2015). Similarly, during the same period of time, the overall rates for gonorrhea and chlamydia in Canada increased by 39 and 58% (75% among men), respectively (see Wilton 2015 and references therein). Interestingly, such a re-emergence follows the dramatic decline in STIs that occurred after the appearance of the HIV in the early 1980s and the consequent widespread use of condoms. However, the introduction of the antiretroviral therapy for HIV in 1996 and a higher adoption of non-condom HIV risk-reduction strategies led to a decreased condom use and the re-emergence of STIs in USA, Canada, and Europe (Fenton and Lowndes 2004; Wilton 2015). In this context, behavioral interventions remain an important tool in the global fight against STIs (Task Force on Community Preventive Services 2007).

For non-reactive SAIS models without decay of awareness, it was shown in Sahneh et al. (2012) that the spread of awareness causes an elevation of the epidemic threshold. Even if the basic reproduction number $R_0 > 1$ in the underlying model of disease transmission, an endemic equilibrium only appears once this elevated threshold is surpassed. However, this result requires the assumption that awareness will not decay over time. In Juher et al. (2015), it was proved that this elevated epidemic threshold disappears if one assumes that awareness will decay over time. Under this assumption, when $R_0 > 1$, awareness may drive the endemic equilibrium to very low levels of disease prevalence, but may not eliminate it or change its stability. Here we show that reactive SAIS models with awareness decay in some cases permit an elevated epidemic threshold so that the endemic equilibrium will disappear and the disease will be driven towards extinction from almost all initial states, even when $R_0 > 1$.

Will a behavioral response that is induced by awareness prevent, all by itself, future flare-ups from low endemic levels? When awareness is assumed to be permanent and demography is not considered, as in Sahneh et al. (2012), Sahneh and Scoglio (2011), the answer is obviously yes. But it is less clear what to expect when awareness decays over time. In the context of ODE models the question translates into one about the existence of sustained oscillations with a significant amplitude. In Sect. 2 we show that such oscillations are ruled out in SAIS models. This result applies to both reactive and non-reactive SAIS models and holds even if we allow more general functional responses rather than rate constants in the models.

However, in Sect. 3 we introduce a more general class of models that we call SAUIS models. They have two distinct compartments of aware hosts who differ in their willingness to alert other susceptible hosts. Thus these models embody the degradation of quality of information as it is transmitted from one individual to another. Our approach for modeling this phenomenon adopts, albeit in greatly simplified form, an idea that was introduced in Agliari et al. (2006) and adopted in Funk et al. (2009) for an epidemic context. We show that sustained oscillations can occur in SAUIS models, even if all rate coefficients are constants. Also, while it seems quite plausible that oscillations could be induced by a time lag between actual prevalence and available information about it, this mechanism is deliberately ruled out by the way we set up our models. Thus our results clearly demonstrate that degradation of information during transmission processes can be the driving mechanism for the existence of periodic waves of infection, supporting the claim in Agliari et al. (2006) that such a degradation can reveal important qualitative and quantitative effects.

The remainder of the paper is organized as follows: In Sect. 2 we formally define reactive SAIS models and prove the results about these models that were described above. In Sect. 3 we formally define reactive SAUIS models and examine some of their basic properties. In particular, we present numerical explorations both of the dynamics in SAUIS models and of regions of the parameter space that are bounded by Hopf bifurcation points. These results reveal intricate possibilities for the dynamics in SAUIS models.

In Sect. 4 we briefly review some related models of behavioral epidemiology that predict oscillations and discuss how they differ from ours. We also discuss some possible implications for public health policy and directions of future work.

2 Reactive SAIS models

2.1 The model

An SAIS model has three compartments: S (susceptible), A (aware) and I (infectious). Susceptible hosts can move to the A-compartment or to the I-compartment, aware hosts can move to the S-compartment due to awareness decay or to the I-compartment due to infection (albeit at a lower rate than susceptible hosts). Upon recovery, infectious hosts can move either to the A-compartment or to the S-compartment.

As was mentioned in the introduction, we will consider awareness as a state of knowledge about the prevalence of the disease that will both induce a behavioral response in the given host and that this host is willing to transmit to other hosts. So, it is natural to assume that awareness can be passed on from aware to unaware individuals like a contagious disease (Epstein et al. 2008). However, while the force of infection in transmission of actual diseases can usually be assumed to be a linear function of the prevalence, the rate at which susceptible individuals become aware usually will show more pronounced nonlinearities. On the one hand, there will be a saturation effect arising from overexposure to information when the proportion of infectious hosts is very high. On the other hand, when the prevalence of the disease is very low, a careless attitude may prevail that renders the channels of transmitting awareness ineffective. Details may significantly vary between different populations and diseases. In order to allow for maximum flexibility in incorporating these effects, we will investigate models where parameters for transmission of awareness are functions of the disease prevalence rather than rate constants. The sources of nonlinearities mentioned above suggest that these parameters may exhibit near switch-like behavior.

The proportions of hosts in the S-, A-, and I- compartments will be denoted by s, a, i respectively. The rates of change of these fractions are governed by the following ODE model:

$$\begin{aligned} \begin{aligned} \frac{da}{dt}&= \alpha _i(i) \, s \, i + \alpha _a(i) \,s \, a + p(i)\, \delta \, i - \beta _a \, a \, i - \delta _a(i) \, a, \\ \frac{di}{dt}&= (\beta \, s + \beta _a\, a - \delta ) \, i, \qquad s+a+i = 1. \end{aligned} \end{aligned}$$

(1)

Here we assume that $\alpha _i(i)$, $\alpha _a(i)$, p(i), $\delta _a(i)$ are nonnegative differentiable functions in [0, 1], $p(i) \le 1$, $\delta _a(0) > 0$, and $\beta , \beta _a, \delta$ are constants such that $0 \le \beta _a < \beta$ and $\delta >0$. Moreover, for $i > 0$ we assume that $\alpha _a(i) > 0$ and $\alpha _i(i) + p(i) > 0$.

The term $\alpha _i(i)i$ represents the rate at which a susceptible host becomes aware due to direct information about the disease prevalence. As $\alpha _i(i)$ may depend on i, the factor i is strictly speaking redundant in this term. But its inclusion simplifies some calculations. Also, inclusion of the factor i makes the similarity with disease transmission more explicit. If $\alpha _i(i)$ is not constant, one can think of the process of obtaining direct information as encountering at least one infectious host and then seeking or re-interpreting independent information about the overall disease prevalence. Thus it seems plausible to assume that $\alpha _i(i)$ is low when infectious hosts are observed so rarely that they are not considered indicative of an outbreak, steeply increases around a critical level of disease prevalence, and then levels off.

Similarly, the term $\alpha _a(i)a$ represents the rate at which susceptible hosts become aware due to a contact with an aware host during which the latter transmits information about the disease. The assumption that $\alpha _a(i) > 0$ is the one that really distinguishes our models from non-reactive SAIS models. The latter can be obtained by simply removing the term $\alpha _a(i) \,s \, a$ from (1). The previously published versions of non-reactive SAIS models also do not include a term $p(i)\, \delta \, i$.

In view of the above discussion we tend to think of $\alpha _i(i)$ and $\alpha _a(i)$ as increasing in a sigmoid-like fashion from 0 or near 0 when $i = 0$ to near a saturation constant when $i = 1$, but these properties are not needed for the results of this section. While one can argue that $\alpha _a(0)$ should be zero so that awareness will not spread in the absence of an outbreak, we also allow for a minimum level of creation of awareness ($\alpha _a(0) > 0$) reflecting, for instance, the spread of rumors and beliefs, or a propensity to become aware because of previous experiences with the disease, even in the absence of current empirical evidence.

The parameter p(i) can be interpreted as the probability that an infectious host will move to the A-compartment as a result of direct experience of the disease. With probability $1-p(i)$ the host will remain oblivious of the dangers posed by the outbreak and will move back to the S-compartment.

The term $\delta _a(i)$ represents the decay of awareness. It could be a constant or any other positive differentiable function, but in realistic models higher disease prevalence should result in slower awareness decay so that $\delta _a(i)$ will be a nonincreasing function of i.

The inequality $\beta _a < \beta$ embodies the assumption that awareness will lead to adoption of a behavioral response that decreases the rate at which hosts contract the infection.

Lemma 1

Assume $\alpha _i(i), \alpha _a(i), p(i), \delta _i(i), \beta$ satisfy the conditions that were spelled out below (1). Then the region $\Omega = \{ (a,i) \in \mathbb {R}^2 \, | \, 0 \le \, a + i \le 1, a \in [0,1] \}$ is forward-invariant.

Proof

By direct inspection of the system we see that $\left. (da/dt)\right| _{a=0} = \alpha _i(i) (1-i) i + p(i)\delta i \ge 0$ for $0 \le i \le 1$, $\left. (di/dt)\right| _{i=0} = 0$, and $\left. (d(a+i)/dt)\right| _{a+i=1} = -(1 - p(i))\delta i - \delta _a(i)a \le 0$. $\square$

2.2 Nullclines and equilibria

By solving $di/dt=0$ we find two parts of the i-nullcline. The first one is given by the horizontal axis $i=0$ and the second one is the straight line

$$\begin{aligned} i(a) = 1 -\frac{\delta }{\beta } - \left( 1 -\frac{\beta _a}{\beta } \right) a, \end{aligned}$$

(2)

which has a slope between $-1$ and 0 under the assumption $\beta _a < \beta$. It intersects the horizontal axis $i=0$ at the point

$$\begin{aligned} a = \frac{\beta - \delta }{\beta - \beta _a}. \end{aligned}$$

By solving $da/dt=0$ we find the a-nullcline. The point (0, 0) always satisfies this equation. For $i > 0$ we have assumed $\alpha _a(i) > 0$ and the other part of this nullcline is implicitly defined by the following equation in the variables i and a:

$$\begin{aligned} a^2 - \left( 1 - i - \frac{(\alpha _i(i) + \beta _a) i + \delta _a(i)}{\alpha _a(i)} \right) a - \frac{\alpha _i(i)}{\alpha _a(i)} \,i \, (1-i) - \frac{p(i)\delta }{\alpha _a(i)}i = 0. \end{aligned}$$

(3)

Since we have assumed $\alpha _i(i) + p(i) > 0$ for $i > 0$, the positive branch of the a-nullcline is given by the graph of the following function a(i) on [0, 1]:

$$\begin{aligned} a(i)= & {} \frac{1}{2} \, \left( 1 - i - \frac{(\alpha _i(i) + \beta _a) i + \delta _a(i)}{\alpha _a(i)} \right. \nonumber \\&\left. + \ \sqrt{\left( 1 - i - \frac{(\alpha _i(i) + \beta _a) i + \delta _a(i)}{\alpha _a(i)} \right) ^2 + 4 \frac{\alpha _i(i)}{\alpha _a(i)} \, i \,(1-i) + \frac{4p(i) \delta }{\alpha _a(i)}i } \,\,\right) . \end{aligned}$$

(4)

Note that $a(1) \ge 0$. If $p(1) = 0$, then $a(1) = 0$. If $p(1) > 0$, the last term under the root in (4) will be positive so that $a(1) > 0$ and the curve a(i) will enter $\Omega$ only at some point (a, i) with $i < 1$. Similarly, $a(0) = 1 - \delta _a(0)/\alpha _a(0)$ if $\delta _a(0) < \alpha _a(0)$, while $a(0) = 0$ if $\delta _a(0) \ge \alpha _a(0)$. Moreover, a(i) is continuous and takes on positive values for all $i \in (0,1)$.

The system has three possible types of equilibria in the first quadrant, namely,

$$\begin{aligned} P_1=(0,0), \quad P_2 = \left( 1-\frac{\delta _a(0)}{\alpha _a(0)}, 0 \right) , \quad P_3 = (a^*, i^*), \end{aligned}$$

where $(a^*, i^*)$ denotes an equilibrium with $i^* > 0$. Since for $i^* > 0$ we have $\alpha _i(i^*) + p(i^*) > 0$ by assumption, any endemic equilibrium must necessarily lie in the interior of $\Omega$. $P_2$ lies in $\Omega$ provided that $\delta _a(0) < \alpha _a(0)$.

Now consider the inequality

$$\begin{aligned} \frac{\beta -\delta }{\beta -\beta _a} > 1 - \frac{\delta _a(0)}{\alpha _a(0)}. \end{aligned}$$

(5)

By sketching the nullclines in the a-i plane, one can see that (5) is a sufficient condition for the existence of at least one endemic equilibrium when $\beta> \delta > \beta _a$, because the function a(i) is continuous and satisfies $a(1) \ge 0$. It will always contain a point on the boundary of $\Omega$ with $a+i = 1$. Thus (5) guarantees that it must intersect the part of the i-nullcline given by (2), which is a straight line such that $i(0) < 1$ with a slope larger than $-1$ and a-intercept $\frac{\beta -\delta }{\beta -\beta _a} < 1$ (see Fig. 1). The condition $\delta > \beta _a$ can be relaxed when all the coefficients in the model are constants; see Lemma 2(a) and its proof.

Fig. 1

Phase portrait of the reactive SAIS with alert rates $\alpha _a(i)=\alpha ^0_a \, (i+1)$ and $\alpha _i(i)=\alpha ^0_i \, (i+1)$, and decay rate $\delta _a(i)=\delta _a^0/(1+i)$ for different values of $\delta _a^0$ showing three possible configurations of equilibria when $R_0 >1$ (top: $\delta _a^0=1$, middle: $\delta _a^0=3$, bottom: $\delta _a^0=5$). Parameters: $p(i) = 0$, $\delta =4$, $\beta =10$, $\beta _a=1$, and $\alpha _a^0=\alpha _i^0=4$. Note that $\alpha _a(0) > 0$ allows the existence of a second equilibrium on the a-axis for small enough values of $\delta _a^0$

Alas, our model assumptions do not rule out multiple endemic equilibria $P_3$. Neither is (5) necessary for the existence of endemic equilibria. Under a number of fairly natural conditions, uniqueness of $P_3$ is guaranteed and (5) is necessary for its existence; see Lemma 2 below. However, as the main results of this section can be derived under our most general assumptions, we will not impose any of these conditions from the outset. We want to mention though that the case of constant rate functions that is covered by points (a1) and (b1) of Lemma 2 is the only one that has been studied in the prior literature on nonreactive SAIS models, and is also the one that directly corresponds to our work in the next section.

Fig. 2

Phase portrait of the reactive SAIS with constant rates showing the existence of two interior equilibria (right) after a saddle-node bifurcation (left) using $\alpha _i$ as a tuning parameter. Parameters: $p=0$, $\beta =6$, $\delta =4$, $\beta _a=2$, $\delta _a=0.9$, $\alpha _a=2$, and $\alpha _i=0.05$ (right) and $\alpha _i=0.1733500838578$ (left)

Lemma 2

Assume $\alpha _i(i), \alpha _a(i), p(i), \delta _i(i), \beta$ satisfy the conditions that were spelled out below (1), and that $\beta > \delta$. Then :

• (a) Under any of the following assumptions, when (5) holds, there exists exactly one interior equilibrium $P_3$:

  • (a1) The functions $\alpha _i(i) = \alpha _i, \alpha _a(i) = \alpha _a, p(i) = p$ and $\delta _a(i) = \delta _a$ are all constant.

  • (a2) The functions $\alpha _i(i), \alpha _a(i), p(i)$ are nondecreasing, $\delta \ge \beta _a$, the function $\delta _a(i)$ is nonincreasing, and $\alpha _a(0) > \beta - \beta _a.$

  • (a3) The functions $\alpha _i(i), \alpha _a(i), p(i)$ are nondecreasing, $\delta \ge \beta _a$, and the function $\delta _a(i)$ decreases steeply enough so that $\beta _a \le - \delta _a'(i)$ for all $i \in [0,1]$.

• (b) If the inequality in (5) is reversed, then any of the following conditions precludes the existence of an interior equilibrium.

  • (b1) The functions $\alpha _i(i) = \alpha _i, \alpha _a(i) = \alpha _a, p(i) = p$ and $\delta _a(i) = \delta _a$ are all constant and

    $$\begin{aligned} \frac{\alpha _a - \delta _a}{\alpha _a + \beta _a} \ge 1 - \frac{\delta }{\beta }. \end{aligned}$$

    (6)

  • (b2) The assumptions of (a2) hold.

  • (b3) The assumptions of (a3) hold.

The proof of this lemma is included in the appendix. In part (b1), without an assumption like (6), saddle-node bifurcations that lead to multiple endemic equilibria are possible; see Fig. 2.

If we define the basic reproduction numbers of the disease and awareness as

$$\begin{aligned} R_0 := \beta /\delta \qquad \text{ and } \qquad R_0^{\,a} := \alpha _a(0)/\delta _a(0), \end{aligned}$$

we can interpret intuitively the conditions for the existence of these equilibria. The disease-and-awareness-free equilibrium $P_1$ always exists. The equilibrium $P_2$ is also disease-free, but has a positive fraction of aware hosts. It exists if, and only if, $R_0^{\,a} > 1$, that is, if in a large and otherwise susceptible population with one aware host, awareness will on average increase. The condition $R_0 > 1$ is necessary and, when $\delta >\beta _a$, condition (5) is sufficient for the existence of an endemic equilibrium $P_3$. Note that condition (5) holds if the disease spreads faster than awareness in the early stages of an outbreak, i.e., if $R_0 > \max \{1, R_0^a\}$. In addition, it also holds when $R_0^a> R_0 > 1$ and $\beta _a$ is close enough to $\delta$, i.e., when the reduction of the transmission probability due to awareness is not significant enough.

The Jacobian matrix of system (1) is

$$\begin{aligned} J = \left( \begin{array}{cc} \alpha _a(i) (s\! -\! a) \! - \! \alpha _i(i) i \! - \!\beta _a i \! -\! \delta _a(i)\ \ &{} \ \ URC\\ - (\beta - \beta _a) i &{} \beta s + \beta _a a - \beta i - \delta \end{array} \right) , \end{aligned}$$

where $URC = (\alpha _i(i)i)'s \!-\! \alpha _i(i)i \!+ \!(\alpha '_a(i) s \!- \!\alpha _a(i) \!- \!\beta _a \!- \!\delta '_a(i)) a \!+ \!(p(i) \delta i)'$ and $s=1-a-i$.

The eigenvalues of J at $P_1$ are

$$\begin{aligned} \lambda _1(P_1) = \alpha _a(0) - \delta _a(0)\qquad \text{ and } \qquad \lambda _2(P_1) = \beta - \delta , \end{aligned}$$

which shows that the necessary conditions for the existence of $P_2$ and $P_3$ imply instability of $P_1$.

At $P_2$, the eigenvalues are

$$\begin{aligned} \lambda _1(P_2) = \delta _a(0) - \alpha _a(0) \qquad \text{ and } \qquad \lambda _2(P_2) = \beta - \delta - (\beta - \beta _a) \left( 1 - \frac{\delta _a(0)}{\alpha _a(0)} \right) . \end{aligned}$$

Note that $\lambda _1(P_2) = - \lambda _1(P_1)$ and that $\lambda _2(P_2) > 0$ if and only if condition (5) holds. The former implies that if $R_0^a > 1$, then $P_1$ is unstable and $P_2$ attracts any trajectory on the a-axis; exactly as one would expect from an SIS-like behavior of awareness.

2.3 Dynamics

Lemma 3

Assume $\alpha _i(i), \alpha _a(i), p(i), \delta _i(i), \beta$ satisfy the conditions that were spelled out below (1). Then the system (1) has no closed orbits inside $\Omega$.

Proof

Let $f_1(a,i)$ and $f_2(a,i)$ denote the functions on the right-hand side of the system. The vector field $\displaystyle (F_1(a,i), F_2(a,i)) = (\frac{1}{a\,i} f_1(a,i), \frac{1}{a\,i} f_2(a,i))$ is $C^1$ in the interior of $\Omega$, and its divergence is given by

$$\begin{aligned} \frac{\partial }{\partial a} F_1(a,i) + \frac{\partial }{\partial i} F_2(a,i) = - \frac{\alpha _i(i)}{a} \left( 1 + \frac{s}{a} \right) - \frac{\alpha _a(i)}{i} - \frac{p(i) \delta }{a^2} - \frac{\beta }{a} < 0 \end{aligned}$$

for all (a, i) in the interior of $\Omega$. So, the divergence does not change sign and does not take the value 0. Therefore, Dulac’s criterion of nonexistence of periodic orbits (Perko 2001) precludes the existence of a closed orbit lying entirely in $\Omega$. $\square$

The next theorem sums up the dynamics of system.

Theorem 1

Assume $\alpha _i(i), \alpha _a(i), p(i), \delta _i(i), \beta$ satisfy the conditions that were spelled out below (1). Then the global behavior of the solutions of the system (1) depends as follows on the remaining parameters:

• (i) If $R_0 \le 1$ and $R^a_0 \le 1$, then $P_1$ is the only equilibrium point and is globally asymptotically stable.

• (ii) If $R_0 \le 1 < R^a_0$, then $P_1$ and $P_2$ are the only equilibrium points. $P_2$ is globally asymptotically stable on $\Omega \backslash \{P_1\}$. When $R_0 < 1$, then $P_1$ is a saddle point.

• (iii) If $R^a_0 \le 1 < R_0$, then no equilibrium $P_2 \ne P_1$ exists in $\Omega$. When $R^a_0 < 1$, then $P_1$ is a saddle point. Each trajectory that starts with $i(0) > 0$ will eventually approach an endemic equilibrium of type $P_3$.

• (iv) If $R_0 > 1$ and $R^a_0 > 1$, then $P_1$ is an unstable point and system (1) has also the equilibrium $P_2$. If $\lambda _2(P_2) < 0$, then $P_2$ is locally asymptotically stable.

• (v) If instead $\lambda _2(P_2) > 0$, then system (1) has also at least one equilibrium $P_3$, with $P_2$ being a saddle point. Each trajectory that starts with $i(0) > 0$ will eventually approach an endemic equilibrium of type $P_3$.

Moreover, when any of the conditions of Lemma 2(a) are satisfied, then the endemic equilibria in point (iii) and (v) are guaranteed to be unique. Similarly, when any of the conditions of Lemma 2(b) are satisfied, then under the assumptions of point (iv) we can conclude that $P_2$ is globally asymptotically stable on $\Omega \backslash \{P_1\}$.

Proof

Recall that $R_0 = \frac{\beta }{\delta }$ and $R^a_0 = \frac{\alpha _a(0)}{\delta _a(0)}$.

For part (i), when $\beta \le \delta$, then the i-nullcline has no points with $s < 1$ and $P_1$ is the only equilibrium point of (1). From Lemmas 1 and 3, together with the Poincaré-Bendixson theorem, it follows that $P_1$ is globally asymptotically stable.

For part (ii), if $\alpha _a(0) > \delta _a(0)$, the a-nullcline intersects the a-axis at $P_2$, and $P_3$ is ruled out as in case (i). So $P_1$ and $P_2$ are the only equilibria of the system. As we have $\lambda _1(P_1) > 0$ and $\lambda _2(P_1) = \beta - \delta \le 0$, it follows that $P_1$ is a saddle point when the inequality $\beta < \delta$ is strict. Moreover, as $(\alpha _i(0) i)' \ge 0$, the first row of the Jacobian at $P_1$ has two nonnegative entries. Thus the eigenvector with eigenvalue $\lambda _1(P_1)$ cannot intersect the interior of $\Omega$. From the same facts that we used in part (i) it follows that $P_2$ is globally asymptotically stable on $\Omega \backslash \{P_1\}$.

Under the assumptions of (iii), we have $\lambda _1(P_1) = \alpha _a(0) - \delta _a(0) \le 0$ and $\lambda _2(P_1) > 0$, and it follows that $P_1$ is a saddle point when the inequality $\alpha _a(0) < \delta _a(0)$ is strict. Moreover, the eigenvector with non-positive eigenvalue lies on the line $i = 0$. Thus at least one endemic equilibrium $P_3$ must exist, and each trajectory that starts with $i(0) > 0$ must eventually approach such an equilibrium.

Under the assumptions of (iv) and (v), the a-nullcline intersects the horizontal axis at $a(0) > 0$ and $P_2$ is always an equilibrium, while $P_1$ is unstable. Recall that the i-nullcline intersects the a-axis at $(\beta - \delta )/(\beta - \beta _a)$ so that the sign of $\lambda _2(P_2)$ determines whether this intersection occurs to the left or to the right of $P_2$. Note that this remains true even when this intersection occurs outside of $\Omega$. The former occurs in case (iv) where $P_2$ is locally asymptotically stable. The latter occurs in case (v), where $P_2$ is a saddle point with the stable direction on the line $i = 0$, while $P_1$ is repelling. It follows that at least one endemic equilibrium $P_3$ must exist, and each trajectory that starts with $i(0) > 0$ must eventually approach such an equilibrium. $\square$

For $R_0^{a} \ge 1$ (so that $P_2 \ge 0$), Theorem 1 shows that the reactive SAIS-model predicts an elevated epidemic threshold which is given by $\lambda _2(P_2)=0$. This threshold can be written by changing the inequality in (5) to an equality. Following the suggestion of one of the reviewers, let us define

$$\begin{aligned} R_0^d = \frac{\beta _a}{\delta } \end{aligned}$$

as the basic reproduction number of aware individuals that would apply to a population consisting entirely of aware hosts. After dividing numerator and denominator of the left-hand side of (5) by $\delta$ and using the definitions of $R_0$ and $R_0^a$, inequality (5) now reads

$$\begin{aligned} \frac{R_0-1}{R_0-R_0^d} > 1 - \frac{1}{R_0^a} = \frac{R_0^a - 1}{R_0^a}. \end{aligned}$$

After passing all terms to the left-hand side of the inequality, multiplying by the positive denominators, canceling $R_0R_0^a$, and then adding 1 to both sides, this can be written equivalently in the usual format for epidemic thresholds as

$$\begin{aligned} R_0 + (R_0^a - 1)(R_0^d-1) > 1. \end{aligned}$$

(7)

When (7) holds and $R_0^a \ge 1$, the disease can always invade the population regardless of how much awareness is initially present, and trajectories that start at an endemic state will eventually approach an endemic equilibrium. In particular, when $R_0^a \ge 1$ and $R_0^d \ge 1$, then (7) holds, as $R_0 = \beta /\delta > R_0^d$ under the assumptions of our model, and the epidemic will spread. On the other hand, if $R_0^a > 1$ but $R_0^d < 1$ (awareness significantly reduces disease transmission), then having $R_0 > 1$ is not enough to guarantee the spread of the epidemic.

When the inequality in (7) is reversed and still $R_0^a \ge 1$, then the disease will not be able to invade a population that has had some prior exposure to awareness, by whatever means, and has reached a state in the basin of attraction of $P_2$. However, it might still be possible for the disease to persist outside of this basin of attraction; see the example of Fig. 2, where the left-hand side of (7) evaluates to 0.888. Therefore, when $R_0^a > 1$ and $R_0^d < 1$, which is feasible in a region of parameter space for our model, (7) defines an elevated epidemic threshold. If in addition any of conditions of Lemma 2(b) hold, then the disease will always be driven to extinction. Under our most general assumptions, the basin of attraction of of $P_2$ will include an open ball around $P_2$ as well as the line segment L that consists of all conditions with $i = 0$ and $a > 0$. It follows that this basin of attraction of $P_2$ must contain an open neighborhood of L. Thus in the limiting case of infinite population size where introduction of a single index case is treated as an initial condition with an infinitesimally small positive i(0), any prior introduction of awareness ($a(0) > 0$) would be sufficient to guarantee that the disease dies out before reaching endemic proportions. This conclusion may not be valid under our general assumptions for small finite populations, but exposure to awareness at a time that is sufficiently early relative to invasion of the disease would still guarantee that it will be driven to extinction under the assumptions of Theorem 1(iv).

It may be of interest to point out that in the limiting case $\beta _a = 0$ the behavioral response that is triggered by awareness confers perfect immunity to becoming infectious. This is quite similar to the assumption that is often made about vaccinations. In this case the inequality in (5) will be reversed if $R_0 < R_0^{a}$. Under these conditions, if all rate functions satisfy the monotonicity conditions of Lemma 2(b3), then any trajectory that starts in $\Omega \backslash \{P_1\}$ will approach the equilibrium $P_2$ with a proportion of $1 - 1/R_0^{a}$ aware hosts, which happens to be exactly the herd immunity threshold for vaccinations that confer perfect protection for diseases with basic reproduction number $R_0^{a}$.

3 SAUIS models

3.1 The model

SAIS models ignore the degradation of information quality as it is transmitted from one individual to another. According to these models, aware individuals would always create aware hosts with the same degree of responsiveness, which seems unrealistic. The following approach to modeling degradation of information was introduced (Agliari et al. 2006) and adopted in Funk et al. (2009) for an epidemic context: We assume that direct information about the disease prevalence induces awareness of the risk of infection and maximum observance of protective measures, while subsequent awareness transmission decreases its impact on an individual’s reaction by a constant decay factor. Following this idea, but simplifying it in order to get a manageable model, we will introduce a new class of individuals whose behavioral response has been induced indirectly through contacts with aware individuals. In contrast to aware individuals, the latter are assumed to be unwilling to convince other people about the risk. They may also show a weaker behavioral response. As they don’t actively participate in the dissemination of awareness, we call them unwilling individuals and let u denote their proportion in the population. Transmission of awareness can create both aware and unwilling hosts.

We will investigate here the simplest equations describing such dynamics that are constructed in direct analogy to the reactive SAIS models of Sect. 2:

$$\begin{aligned} \begin{aligned} \frac{da}{dt}&= \alpha _i \, s \, i + \alpha _a \, s \, a + p \,\delta \, i - \beta _a \, a \, i - \delta _a\, a, \\ \frac{du}{dt}&= \delta _a\, a + \alpha _u \, s \, a + q \, \delta \, i - \beta _u\, u \, i - \delta _u\, u, \\ \frac{di}{dt}&= (\beta \, s + \beta _a \, a + \beta _u \, u - \delta )\, i, \qquad s+a+u+i = 1. \end{aligned} \end{aligned}$$

(8)

Here $\alpha _u\, a$ is the rate at which susceptible hosts become unwilling after having a contact with an aware host, and $\delta _u$ is the rate of awareness decay of the unwilling hosts. This model implicitly assumes that aware hosts first turn into unwilling hosts before possibly entering the susceptible compartment. It also allows for coexistence of three behavior patterns at the time of recovery from the disease: The host will then move into the A-compartment with probability p, into the U-compartment with probability q, and into the S-compartment with probability $1-p-q$. All other terms play the same role as the corresponding terms in the reactive SAIS model.

We assume that all rate coefficients are constants and, with the possible exception of $\beta _a$ and $\beta _u$, are positive, with $0 \le \beta _a, \beta _u < \beta$. Similarly to the SAIS model, one could allow some of the rate coefficients to depend on i; see Sect. 4 for a brief discussion of these extensions of the model. However, we specifically restrict our attention here to the case of constant rate coefficients to emphasize the point that their dependence on i is not a necessary condition for observing periodic oscillations.

Lemma 4

The region $\Omega = \{ (a,u,i) \in \mathbb {R}^3_+ \, | \, 0 \le \, a + u + i \le 1\}$ is positively invariant.

Proof

Direct inspection of the system (8) shows that $\left. (da/dt)\right| _{a=0} \ge 0$ and the inequality is strict when $si > 0$. Similarly, $\left. (du/dt)\right| _{u=0} \ge 0$, and the inequality is strict when $a > 0$. On the other hand, $\left. (di/dt)\right| _{i=0} = 0$. It follows that the (a, i)- and (u, i)-coordinate planes repel the trajectories and that the (a, u)-plane is invariant. So, we only need to see that trajectories cannot cross the boundary $a+u+i=1$. If $\mathbf {v}$ denotes the vector field defined by the right-hand side of the system, we can see this by computing the scalar product of $\mathbf {v}$ on this boundary with the outward-pointing normal vector $\mathbf {n}=(1,1,1)$. Since $s = 0$ in this region of the boundary, we get $\mathbf {v} \cdot \mathbf {n} = - \delta _u \, u - \delta \, i \, (1 - p - q) \le 0$. The inequality is strict except at the point (1, 0, 0), where $da/dt = - du/dt < 0 = di/dt$. Therefore the vector field on the boundary $a+u+i=1$ never points towards the exterior of $\Omega$. $\square$

3.2 Possible equilibria

The system (8) can have up to three types of equilibria in $\Omega$.

The first one is the disease-free, awareness-free equilibrium $P_1 = (0,0,0)$.

The second one is the disease-free equilibrium $P_2 = (a^*_0, u^*_0, 0)$. From the first line of (8) we have $s^*_0 = \delta _a/\alpha _a$. In view of the equality $u^*_0 = 1 - a^*_0 - \delta _a/\alpha _a$, the second line of (8) implies:

$$\begin{aligned} a^*_0 = \frac{ \delta _u \left( 1 - \frac{\delta _a}{\alpha _a} \right) }{ \delta _a \left( 1 + \frac{\alpha _u}{\alpha _a} \right) + \delta _u}, \quad u^*_0 = \left( 1 - \frac{\delta _a}{\alpha _a} \right) \frac{\delta _a \left( 1 + \frac{\alpha _u}{\alpha _a} \right) }{\delta _a \left( 1 + \frac{\alpha _u}{\alpha _a} \right) + \delta _u}. \end{aligned}$$

(9)

The third kind of possible equilibrium is an endemic equilibrium, i.e., a point $P_3=(a^*,u^*,i^*)$ of $\Omega$ with

$$\begin{aligned} i^* = 1 - \left( 1-\frac{\beta _a}{\beta } \right) a^* - \left( 1 - \frac{\beta _u}{\beta } \right) u^* - \frac{\delta }{\beta } > 0. \end{aligned}$$

(10)

As the upper panel of Fig. 3 shows, there may be more than one endemic equilibrium in the interior of $\Omega$.

3.3 Existence and linear stability of equilibria

The disease-free, awareness-free equilibrium $P_1 = (0,0,0)$ always exists. Evaluating the Jacobian matrix of system (8) at $P_1$ we have

$$\begin{aligned} J(P_1) = \left( \begin{array}{ccc} \alpha _a - \delta _a \ \ &{} \ \ 0 \ \ &{} \ \ \alpha _i + p \delta \\ \delta _a + \alpha _u \ \ &{} \ \ - \delta _u \ \ &{} \ \ q \delta \\ 0 \ \ &{} \ \ 0 \ \ &{} \ \ \beta - \delta \end{array} \right) , \end{aligned}$$

(11)

whose eigenvalues are

$$\begin{aligned} \lambda _1(P_1)=\alpha _a - \delta _a, \quad \lambda _2(P_1)= - \delta _u, \quad \lambda _3( P_1)= \beta - \delta . \end{aligned}$$

So, as expected, when $\beta > \delta$, i.e., when $R_0 > 1$, this equilibrium is unstable. Moreover, as in the reactive SAIS model, when $R^a_0 := \alpha _a/\delta _a > 1$, this equilibrium is unstable independently of the sign of $\beta - \delta$ and $P_2$ becomes biologically meaningful.

By (9), $P_2$ exists in $\Omega \backslash \{P_1\}$ if, and only if, $R^a_0 = \alpha _a/\delta _a > 1$. In this case it must be in the interior of the intersection of the a-u-plane with $\Omega$. The Jacobian matrix of (8) at this equilibrium is

$$\begin{aligned} J(P_2) = \left( \begin{array}{ccc} - \alpha _a\, a^*_0 \ \, &{} \ \, -\alpha _a\, a^*_0 \ \, &{} \ \, \alpha _i s^*_0 - (\alpha _a+\beta _a)a^*_0 + p \delta \\ \delta _a + \alpha _u (s^*_0 - a^*_0) \ \, &{} \ \, -\alpha _u a^*_0 - \delta _u \ \, &{} \ \, -\alpha _u a^*_0 - \beta _u u^*_0 + q \delta \\ 0 \ \, &{} \ \, 0 \ \, &{} \ \, \beta s^*_0 + \beta _a a^*_0 + \beta _u u^*_0 - \delta \end{array} \right) ,\qquad \end{aligned}$$

(12)

with $s^*_0 = \delta _a/\alpha _a$. Here we have used the observation that the expression $\alpha _a\, (s_0^* - a^*_0) - \delta _a$ obtained from direct computation of the partial derivative in the upper left corner of $J(P_2)$ simplifies to $-\alpha _a a^*_0$. When one computes the determinant of the submatrix formed by the intersection of the first two rows of $J(P_2)$ with its first two columns, the only negative term cancels out. Moreover, the trace of this submatrix is negative, so that $J(P_2)$ has two eigenvalues that are either negative or have negative real parts.

The third eigenvalue $\lambda _3(P_2)=\beta - (\beta - \beta _a) \, a^*_0 - (\beta - \beta _u) \, u^*_0 - \delta$ is negative if $\beta < \delta$ and $P_2 \in \Omega$. When $\beta > \delta$ and

$$\begin{aligned} \frac{\beta - \beta _a}{\beta - \delta } \, a^*_0 + \frac{\beta - \beta _u}{\beta - \delta } \, u^*_0 < 1, \end{aligned}$$

(13)

then $\lambda _3(P_2)$ will be positive. By (9), the values of $a_0^*, u_0^*$ do not depend on the disease transmission parameters, and it can be seen from (13) that there are large regions of the parameter space where $\lambda _3(P_2)$ is positive and large regions where it is negative while $\beta > \delta$.

3.4 Transcritical bifurcations

3.4.1 Transcritical bifurcation at $R^a_0=1$

Assume $\beta < \delta$ so that $\lambda _3(P_1) < 0$. As $\lambda _1(P_1)$ changes from negative to positive when the bifurcation parameter $R^a_0=\frac{\alpha _a}{\delta _a}$ increases past 1, the equilibrium $P_1$ loses its stability at the bifurcation value 1. Simultaneously, $P_2$ enters the biologically feasible region $\Omega$ and becomes locally asymptotically stable as explained in the previous subsection. Moreover, at the bifurcation value $P_1$ and $P_2$ coincide. Note also that the plane $i = 0$ is invariant for the system (8), and when $a^*_0$ is negative, then the determinant of the submatrix formed by the intersection of the first two rows of $J(P_2)$ with its first two columns is negative. This implies that when $P_2$ crosses into the biologically feasible region, the equilibria $P_1$ and $P_2$ interchange their stability.

3.4.2 Transcritical bifurcation at $R_0=1$

For the analysis of the bifurcations of the endemic equilibrium $P_3$ from $P_1$ and $P_2$, we will use standard results for the existence of a transcritical bifurcation [see the criterion that is given right after Sotomayor’s Theorem in Perko (2001)]. In order to use them, we introduce some notation. Let $\mathbf f$ denote the vector defined by the right-hand side of system (8) and let $\mathbf{f}_{\mu }$ be the vector of partial derivatives of its components $f_i$ with respect to a bifurcation parameter $\mu$. Moreover, let $D\mathbf{f}_{\mu }$ be the Jacobian matrix of $\mathbf{f}_{\mu }$ and let $D^2\mathbf{f}(\mathbf{y}, \mathbf{y})$ be the column vector with components $\left( D^2\mathbf{f}(\mathbf{y}, \mathbf{y}) \right) _k := \sum _{j,l} \frac{\partial ^2 f_k}{\partial x_j \partial x_l} y_j y_l$, where $\mathbf y$ is a vector in $\mathbb {R}^3$, $x_1=a$, $x_2=u$, and $x_3=i$. We will use $\mathbf{f}^{BP}, \mathbf{f}^{BP}_\mu$ to indicate that the relevant objects are computed at the bifurcation point.

The endemic equilibrium $P_3$ can bifurcate from $P_1$ when $\alpha _a < \delta _a$ (i.e., $R_0^a <1$). In particular, since $\lambda _3(P_1)=\beta -\delta$ is a simple eigenvalue, a bifurcation occurs for $\beta = \delta$ (i.e., at $R_0 =1$). Taking $\beta$ as the bifurcation parameter and evaluating the Jacobian matrix $J(P_1)$ at the bifurcation point, it follows that the row vector $\mathbf{u} = (0,0,1)$ and the column vector $\mathbf{v} = (1, (\delta _a + \alpha _u + \frac{q\delta (\delta _a - \alpha _a)}{\alpha _i + p\delta })/\delta _u, (\delta _a - \alpha _a)/(\alpha _i + p\delta ))^T$ are the left and right eigenvectors for $\lambda _3=0$, respectively. Moreover, $\mathbf{f}_\beta = (0,0, (1-a-u-i)i)^T$. A straightforward computation at the bifurcation point leads to

• (1) $\mathbf{u} \cdot \mathbf{f}_{\beta }^{BP} = 0$,

• (2) $\mathbf{u} \cdot (D\mathbf{f}_{\beta }^{BP} \mathbf{v}) = v_3 = (\delta _a - \alpha _a)/(\alpha _i + p \delta ) > 0$, and

• (3) $\mathbf{u} \cdot \left( D^2 \mathbf{f}^{BP}(\mathbf{v}, \mathbf{v}) \right) = -2v_3 \left( (\beta -\beta _a) v_1 + (\beta -\beta _u) v_2 + \beta v_3 \right) < 0$.

Note that the inequality $\alpha _a < \delta _a$ and the assumption $0 \le \beta _a, \beta _u < \beta$ of the SAUIS model imply that all coordinates of $\mathbf v$ and coefficients in (3) are positive, so that the whole expression becomes negative and, in particular, nonzero. So, when $R^a_0 < 1$ system (8) experiences a transcritical bifurcation as $\beta$ crosses the bifurcation value $\beta =\delta$ (Perko 2001). Moreover, Theorem 4.1 in Castillo and Song (2004), together with the inequality $> 0$ in (2) and the inequality $< 0$ in (3), implies that the direction of the bifurcation is always the same, namely, system (3) experiences a forward bifurcation at $R_0 = 1$.

3.4.3 Transcritical bifurcation at $\lambda _3(P_2)=0$

Let us now assume $\alpha _a > \delta _a$ to guarantee the existence of a positive $P_2$, and $\beta > \delta$ to allow $\lambda _3(P_2)$ to be positive for some parameters values. From the discussion surrounding (13) it follows that $\lambda _3(P_2)=0$ if and only if

$$\begin{aligned} \frac{\beta - \beta _a}{\beta - \delta } \, a^*_0 + \frac{\beta - \beta _u}{\beta - \delta } \, u^*_0 = 1, \end{aligned}$$

(14)

with $a^*_0$ and $u^*_0$ given by (9). That is, at this parameter combination, the endemic equilibrium $P_3$ bifurcates from $P_2$, the disease-free equilibrium with aware individuals.

When the Jacobian matrix (12) is evaluated at the bifurcation point, the value in the lower right corner becomes zero. We get

$$\begin{aligned} J^{BP} = \left( \begin{array}{ccc} - \alpha _a\, a^*_0 \ \, &{} \ \, -\alpha _a\, a^*_0 \ \, &{} \ \, \alpha _i \delta _a/\alpha _a - (\alpha _a+\beta _a)a^*_0 + p \delta \\ \delta _a + \alpha _u (\delta _a/\alpha _a - a^*_0) \ \, &{} \ \, -\alpha _u a^*_0 - \delta _u \ \, &{} \ \, -\alpha _u a^*_0 - \beta _u u^*_0 + q \delta \\ 0 \ \, &{} \ \, 0 \ \, &{} \ \, 0 \end{array} \right) , \end{aligned}$$

which has the row vector $\mathbf{u} = (0,0,1)$ as left eigenvector with eigenvalue $\lambda _3 = 0$. Since the determinant $\det \left( J_2^{BP} \right) = \alpha _a a^*_0(\delta _a+\delta _u+\alpha _u\delta _a/\alpha _a) > 0$ of the submatrix $J_2^{BP}$ of $J^{BP}$ formed by the first two rows and the first two columns is strictly positive, the first two columns of $J_2^{BP}$ are linearly independent. Therefore, the third component $v_3$ of the right eigenvector $\mathbf v$ with eigenvalue $\lambda _3=0$ must be different from 0.

If we take $\mu =\beta _a$ as the bifurcation parameter, then $\mathbf{f}_{\beta _a} = (-a\,i, 0, a\,i)^T$. As before, a straightforward computation at the bifurcation point leads to:

• (1) $\mathbf{u} \cdot \mathbf{f}_{\beta _a}^{BP} = 0$,

• (2) $\mathbf{u} \cdot (D\mathbf{f}_{\beta _a}^{BP} \mathbf{v}) = a^*_0 v_3 \ne 0$, and

• (3) $\mathbf{u} \cdot \left( D^2 \mathbf{f}^{BP}(\mathbf{v}, \mathbf{v}) \right) = -2v_3 \left( (\beta -\beta _a) v_1 + (\beta -\beta _u) v_2 + \beta v_3 \right)$.

Recall that $\mathbf v$ is by definition orthogonal to the first two rows of $J^{BP}$ and these are linearly independent. Thus $\mathbf v$ cannot also be orthogonal to any vector that is linearly independent of the latter two vectors. Therefore, since $v_3 \ne 0$, the inequality $\mathbf{u} \cdot \left( D^2 \mathbf{f}^{BP}(\mathbf{v}, \mathbf{v}) \right) \ne 0$ will follow whenever the vector given by $\mathbf{b} := (\beta -\beta _a, \beta -\beta _u, \beta )$ is linearly independent of the two first rows of $J^{BP}$. This may not always be the case, but it will be generically true. For instance, the first two columns of $J^{BP}$ are linearly independent and the parameter $\alpha _i$ appears in only one position of $J^{BP}$, while $a_0^*$ and $u_0^*$ that give the location of $P_2$ do not depend on it. Since it does not appear in $\mathbf b$, one can take $\alpha _i$ as a free parameter in order to see whether the previous inequality fails for a particular positive choice of $\alpha _i$ given any choice of all other parameters.

Thus the criterion given in Perko (2001) after Sotomayor’s Theorem guarantees that for a nonempty open set of parameter settings at which $P_1 \ne P_2 \in \Omega$ the system (8) experiences a transcritical bifurcation as $\beta _a$ passes through the bifurcation value

$$\begin{aligned} \beta _a^c :=\beta - \frac{1}{a_0^*} \left( \beta -\delta - (\beta -\beta _u)u_0^* \right) . \end{aligned}$$

However, in contrast to what happens at $R_0=1$, the direction of the bifurcation is not always the same. To illustrate this fact, Fig. 3 shows an example of forward and backward bifurcations occurring at $\lambda _3(P_2)=0$ for $\beta _a < \beta$. Here $p = q = 0$, but similar examples with $p, q > 0$ exist.

Fig. 3

Transcritical bifurcation diagrams of system (8) for $\beta =2$, $\beta _u=1$, $\delta =1$, $\delta _a=0.01$, $\delta _u=0.05$, $\alpha _i=0.8$, $\alpha _u=0.1$, $p = q = 0$, $\alpha _a=0.1$ (top) and $\alpha _a=1$ (bottom). The stable (unstable) equilibria are depicted with a solid (dashed) line. Bifurcation values: $\beta _a^c=0.844$ and $\beta =0.62$ (for the fold bifurcation) (upper panel); $\beta _a^c=0.988$ (lower panel)

The same conclusion holds if we use $\beta$ or $\beta _u$ as a bifurcation parameter. However, the use of $\beta _a$ seems more appropriate because its value is related to the effectiveness of the response adopted by aware hosts.

3.5 Hopf bifurcations

We will see that in contrast to the SAIS model, sustained oscillations are possible in the SAUIS model thanks to the occurrence of Hopf bifurcations. A pair $(\sigma , \tau )$ of parameters of (8) will be called a Hopf pair if there exists an equilibrium point at which the Jacobian matrix has a pair of pure imaginary eigenvalues. This definition of course depends on chosen values of the other parameters. For simplicity we will always implicitly assume that these other parameters are fixed and suppress them in our notation.

An explicit criterion that specifies whether an $n \times n$ matrix M, with coefficients that may depend upon parameters, has a pair of pure imaginary eigenvalues is given in Guckenheimer et al. (1997). For our purposes the case when $n = 3$ and $M=J$ is the Jacobian matrix at endemic equilibrium $P_3$ is relevant. Let $p(\lambda )=\lambda ^3 + c_2 \lambda ^2 + c_1 \lambda + c_0$ be the characteristic polynomial of J. Its coefficients are $c_0 = -\det (J)$, the sum  $c_1$ of the principal minors of J, and $c_2 = - \text {trace}(J)$. According to Theorem 2.1 and Table 1 in Guckenheimer et al. (1997), the matrix J has precisely one pair of pure imaginary eigenvalues if and only if

$$\begin{aligned} c_0 - c_1 c_2 = 0 \quad \text {and} \quad c_1 > 0. \end{aligned}$$

(15)

Therefore, to locate the Hopf pairs in a $(\sigma , \tau )$ parameter space, we will use $\sigma$ as free parameter and solve the system given by the equilibrium equations

$$\begin{aligned}&\alpha _i (1-a^*-u^*-i^*) i^* + \alpha _a (1-a^*-u^*-i^*) a^* - \beta _a a^* i^* - \delta _a a^* + p \delta i^* = 0 \nonumber \\&\quad \delta _ a a^* + \alpha _u (1-a^*-u^*-i^*) a^* - \beta _u u^* i^* - \delta _u u^* + q \delta i^* = 0, \end{aligned}$$

(16)

combined with (10) and (15), for $a^*$, $u^*$, $i^*$, and $\tau$, while all other parameters are set to fixed values.

The component $\tau$ of the solution, if any, and the corresponding value of $\sigma$ define a Hopf pair $(\sigma , \tau )$ for which a Hopf bifurcation occurs at the endemic equilibrium $P_3=(a^*,u^*,i^*)$. The set of Hopf pairs defines the so-called Hopf-bifurcation curve H in $(\sigma , \tau )$ parameter space,

$$\begin{aligned} H = \{ (\sigma , \tau ) \in \mathbb {R}^2_+ \, | \, \exists P_3 \in \mathbb {R}^3_+ \, \text {for which} \, ~(10), \, (15), \, \text {and} \, (16) \, \text {hold} \}. \end{aligned}$$

Note that H can be parametrized by $i^*$ (or any of the components of $P_3$) (Szabó et al. 2012).

We will focus here on the case $\sigma = \beta _a$ and $\tau = \alpha _i$. This choice appears to be of particular interest, as $\alpha _i$ may be most amenable to alteration by increased epidemiological monitoring (see Sect. 4), while $\beta _a$ can be thought of as inversely proportional to the effectiveness of the behavioral response in aware hosts (see Sect. 4). An example of a curve H in the $(\beta _a, \alpha _i)$ parameter space is shown in Fig. 4. For each value of $\beta _a$ within the range [0, 0.4379), there exist two solutions of Eqs. (15) and (16). The prevalence of the disease at the bifurcation points along the curve H is presented in Fig. 5. The figure clearly shows how $i^*$ decreases monotonously with the value of $\alpha _i$ at the Hopf pairs.

Fig. 4

Hopf-bifurcation curve H of system (8) for $\delta =1$, $\delta _a=0.01$, $\delta _u=0.05$, $\beta =3$, $\beta _u=0.5$, $\alpha _a=0.01$, $\alpha _u=1$, and $p = q = 0$. For pairs $(\beta _a,\alpha _i)$ inside the region bounded by this curve and the $\alpha _i$-axis system (8) has an unstable endemic equilibrium and a stable periodic orbit

The (sub- or supercritical) character of the Hopf bifurcation can be found by computing the first Lyapunov exponent, which is close to $-2.5$ for intersection of the straight line $\beta _a=0.2$ with the upper branch of the curve of predicted bifurcation points and close to $-0.25$ for the intersection of this line with the lower branch. Thus for each pair $(\beta _a,\alpha _i)$ inside $\mathcal {R}$, the endemic equilibrium $P_3$ is unstable, as $J(P_3)$ has two complex eigenvalues with positive real part and a third eigenvalue that is real and negative, while stable periodic orbits appear around $P_3$. Outside this region, $P_3$ is asymptotically stable. Therefore, the points on the curve H define supercritical Hopf bifurcations of the system (8). The numerical explorations in Fig. 6 confirm this prediction.

Fig. 5

Fraction of infectious hosts as a function of $\alpha _i$ along the Hopf-bifurcation curve H of system (8) for $\delta =1$, $\delta _a=0.01$, $\delta _u=0.05$, $\beta =3$, $\beta _u=0.5$, $\alpha _a=0.01$, $\alpha _u=1$, and $p = q = 0$. At both ends of the curve, $\beta _a=0$

Numerical simulations confirm our predictions regarding the Hopf pairs shown in Fig. 4. In particular, for $\beta _a = 0.2$ and the parameter settings as specified in the caption of Fig. 4, Hopf bifurcations in system (8) are predicted for $\alpha _i^* = 0.2251$ and $\alpha _i^{**} = 0.8916$. Figure 6 illustrates what happens if we increase $\alpha _i$ from below $\alpha _i^*$ to above $\alpha _i^{**}$ for trajectories with initial condition $a(0)=u(0)=0$ and $i(0)=0.1$, far from the endemic equilibrium $P_3$.

The top left panel of Fig. 6 shows a simulation with $\alpha _i = 0.2$, right below the first Hopf bifurcation. The solution quickly approaches an endemic equilibrium with disease prevalence $i^* = 0.1249$. When we increase $\alpha _i$ slightly above $\alpha _i^*$, to 0.24, then we observe significant oscillations in the top right panel. Moreover, these oscillations correspond to a stable limit cycle in the phase portrait of the system that attracts other trajectories whose initial conditions are not necessarily close to the endemic equilibrium. For the (unstable) endemic equilibrium we get $i^* = 0.0651$ in this case. Due to the oscillations, the long-term mean prevalence will be even slightly lower, around 0.0592.

A similar picture of persistent oscillations is shown in the bottom left panel of Fig. 6 where we chose $\alpha _i = 0.5$, about half-way between $\alpha _i^*$ and $\alpha _i^{**}$. The existence of a stable limit cycle becomes clearly noticeable in simulations over a longer time horizon. Here the prevalence at endemic equilibrium is $i^* = 0.0148$, very close to the mean prevalence $\overline{i} = 0.0150$ in the long run. Note that for these settings long periods of very low disease prevalence alternate with short spikes that indicate periodic small flare-ups. These flare-ups are followed by a rapid increases in awareness, which indicate a panic-like spread of information.

Finally, when we further increase $\alpha _i$ beyond $\alpha _i^{**}$ to 0.94, we initially observe a similar pattern of alternating periods of extremely low disease prevalence, interrupted by small flare-ups with panic-like spread of awareness (see the bottom right panel in Fig. 6). The prevalence at endemic equilibrium further decreases to $i^* = 0.0065$. However, as this panel indicates, the amplitude of these oscillations now decreases, albeit very slowly, because the pair $(\beta _a, \alpha _i)=(0.2, 0.94)$ is very close to the curve H.

Fig. 6

Evolution of the fraction of infectious (solid line), aware (dashed line), and unwilling (dot-dashed line) hosts according to system (8) for different values of $\alpha _i$ along the vertical section in Fig. 4 corresponding to $\beta _a=0.2$: $\alpha _i=0.2$ (top left), 0.24 (top right), 0.5 (bottom left), 0.94 (bottom right). Fixed parameters: $\alpha _a = 0.01, \ \alpha _u = 1, \ \delta = 1, \ \delta _a = 0.01, \ \delta _u = 0.05, \ \beta = 3, \ \beta _a = 0.2, \ \beta _u = 0.5, \ p = q = 0$. Initial condition: $a(0) = u(0) = 0$, $i(0) = 0.1$

It is interesting to compare Fig. 6 with Fig. 7. The parameter settings in the latter are similar to those in the upper left panel of Fig. 6, except that now p and q assume small positive values. In particular, compare the upper left panels of these figures. While sustained oscillations are absent when $p = q = 0$ (Fig. 6), they do occur when $p = 0.05$ and q is very small. As the two lower panels of Fig. 7 show, increasing q first dampens and then eliminates these oscillations. Thus we conclude that for these settings of the remaining parameters the oscillations are driven by having a positive proportion of hosts who move into the A-compartment upon recovery. The upper right panel of Fig. 7 corresponds to the upper right panel of Fig. 6 and shows similar oscillations, but with decreased amplitude when $p, q > 0$.

While p and q are fairly small in the parameter settings for Fig. 7, we found that sustained oscillations are possible even when $p+q = 1$. Figure 8 shows the Hopf bifurcation curve for the Hopf pair (p, q) for a similar setting of the other parameters. The first Lyapunov exponent for the point on the Hopf bifurcation curve with $p = 0.1$ is close to $-2.6$. Thus sustained oscillations occur in the area under the curve. The line segment inside this region indicates locations where $p+ q = 1$. Interestingly enough, for q around 0.18, the curve predicts no oscillations for values of p that are very small or very close to 1, while it predicts sustained oscillation when p takes moderate values. For more explorations of the influence of parameters p and q on the dynamics of the model, see Xin (2016).

Fig. 7

Evolution of the fraction of infectious (solid line), aware (dashed line), and unwilling (dot-dashed line) hosts for different values of q corresponding to $p=0.05$; $q=0.05$ (top left), 0.1 (bottom left), 0.2 (bottom right). Initial condition: $a(0) = u(0) = 0$, $i(0) = 0.1$. In the left and bottom right panels, $\alpha _u = 1$, $\delta = 1$, $\delta _a = 0.01$, $\delta _u = 0.05$, $\beta = 3$, $\beta _u = 0.5$, $\alpha _a = 0.01$, $\beta _a=0.2$ and $\alpha _i=0.2$. In the top right panel, $\alpha _i = 0.24$, $p = 0.1$, $q = 0.1$, and all the other parameters stay the same

Fig. 8

Hopf-bifurcation curve H of system (8). The straight line sets the boundary $p + q = 1$, so only the part of the Hopf-bifurcation curve that lies on the left of this line makes biological sense. When $p = 0.05$, $q^*=0.1923$ is the only Hopf-bifurcation point. Here, $\alpha _u = 3$, $\delta = 1$, $\delta _a = 0.01$, $\delta _u = 0.05$, $\beta = 3$, $\beta _u = 0.5$, $\alpha _a = 0.012$, $\beta _a = 0.2$, $\alpha _i = 0.05$

4 Discussion

Previous results on various types of SAIS models had suggested that a behavioral response of hosts who are aware of an ongoing outbreak of a disease can be effective in eliminating the endemic equilibrium or driving it to very low levels. For non-reactive SAIS models ($\alpha _a = 0$) it was shown in Sahneh et al. (2012) that if there is no awareness decay $(\delta _a = 0)$, then there is an elevated epidemic threshold so that all outbreaks will be minor as long as $\beta _a< \delta < \beta$ and $\alpha _i$ is sufficiently large. In Juher et al. (2015), it was proved that this elevation of the epidemic threshold disappears if one assumes that awareness will decay over time. The reason is that in order to suppress epidemic spreading, we need a permanent fraction of aware hosts. If awareness is eventually lost, then the presence of aware individuals created at the early stages of an epidemic is not enough to contain an outbreak, and the classic disease-invasion epidemic threshold $R_0 = \beta /\delta = 1$ drives the dynamics.

Here we have introduced and studied reactive SAIS models that make the natural assumption that awareness can also be transmitted from an aware host to a susceptible one ($\alpha _a > 0$). We have shown that under certain conditions in these models a permanent fraction of aware individuals can be sustained even in the presence of awareness decay, which again can lead to elevation of the epidemic threshold above $R_0$.

The question of when the response will be sufficient to prevent future flare-ups from low endemic levels had not previously been addressed in the literature. This question is of interest for models where awareness will decay over time. Lemma 3 of Sect. 2 indicates that when there is not much differentiation between hosts in their propensity to share information about the disease with other hosts, future flare-ups are ruled out. This result applies both to SAIS models with constant rates and to SAIS models in which the rate coefficients depend on the prevalence of the disease.

However, in most real human populations hosts will significantly differ in how effectively they contribute to growing awareness of other hosts. We constructed a new class of models, SAUIS models, that incorporate this phenomenon. We think of them as the simplest possible straightforward generalization of SAIS models in this direction. SAUIS models permit persistent cycles of flare-ups that induce panic-like spread of information, which will temporarily drive the prevalence to low levels without leading to elimination of the disease or providing permanent protection against future major outbreaks. What we observe in such dynamics are flare-ups of the disease, closely followed by steep rises of hosts who are willing to spread awareness, which in turn are followed by increases in the number of unwilling hosts and corresponding decreases in the size of the S-compartment. This drives the disease prevalence to low levels while the numbers of unwilling hosts stay high, but the spread of awareness slows down until the S-compartment gets sufficiently replenished for the next flare-up of the disease.

This phenomenon does not depend on variable rate coefficients, as all of these are assumed constant in the version of SAUIS models presented here. It can occur regardless of whether some or all hosts are assumed to move into the union of the A-compartment and the U-compartment at the time of recovery from the disease. While increasing the fraction of hosts who move to the U-compartment tends to dampen and eventually prevent oscillations, increasing the fraction of hosts who move into the A-compartment as a result of direct experience can sometimes induce and sometimes prevent oscillations. This, together with the contrasting result for the closely related SAIS models, clearly demonstrates that the observed oscillations are in fact driven by the unequal propensity of hosts to share information when differences in the willingness to share information arise from the degradation of the information during subsequent transmission events as modeled in Funk et al. (2009).

These findings seem obviously relevant from the point of view of designing effective policies for controlling infections. As the two top panels of Fig. 6 suggest and our calculations confirmed, the average disease load, as interpreted as the mean prevalence calculated from the area under the solid curve, may in the long run be similar and sometimes even lower in the presence of periodic flare-ups than for similar parameter settings where an endemic equilibrium is approached. Thus if the primary goal of control is a reduction of average load, then elimination of periodic flare-ups may sometimes even be counterproductive. However, peak load may be more important than average load in terms of the danger of overwhelming the health care system (Fefferman 2014), and it will be higher under oscillatory dynamics. Thus if the primary goal of control is to reduce peak load of the disease as much as possible, then control measures should be targeted at elimination of possible flare-ups.

One possible control strategy is to commit resources to continuous monitoring the prevalence i of the disease and dissemination of this information by health care professionals and the media. In terms of SAUIS models, this can be interpreted as increasing the parameter $\alpha _i$. It seems plausible that this strategy will decrease the value of $i^*$ at the endemic equilibrium and the long-term average disease load $\bar{i}$. Our numerical explorations suggest that this intuition is correct. When $\alpha _i$ crosses the threshold represented by the lower branch of the curve in Fig. 4 the endemic equilibrium will become unstable, with oscillatory dynamics for values above the threshold. However, as Fig. 6 shows, right above the threshold the resulting oscillations may result in an increased peak load relative to values of $\alpha _i$ right below it. A further increase of $\alpha _i$ beyond a second threshold represented by the upper branch of the curve in Fig. 4 will lead to renewed local asymptotic stability of the endemic equilibrium. But as the lower-right panel Fig. 6 shows, even above this threshold damped oscillations may remain observable for a considerable time horizon when trajectories start far away from the equilibrium. These and other results presented in Sect. 3 reveal a surprisingly rich and intricate dynamics of SAUIS models. It becomes clear that optimizing the dynamics by controlling one or more parameters that govern the spread of awareness will usually be quite delicate even if one could treat the model at face value.

Our findings open several avenues of future research. It would be interesting to analytically derive necessary and/or sufficient conditions on the information flow about the disease between different types of hosts that would preclude the existence of sustained oscillations. However, in order to base actual public health policy on such results, they would need to be sufficiently robust and carry over to wider classes of models of transmission of awareness and degradation of information. The SAUIS models that we introduced and explored here are useful for clearly demonstrating effects that are precluded by the simplifying assumptions of previously studied SAIS models. However, in biologically more realistic models one would presumably want to consider more fine-grained scenarios of degradation of information, perhaps by incorporating a chain of compartments $A_1, A_2, \dots , A_m$ that would represent progressively decreasing willingness to share the information and/or decreasing strength of the behavioral response. Such chains could also be used to more directly adapt the model of Funk et al. (2009).

Moreover, the assumption of constant rate coefficients is clearly an oversimplification. We already discussed in the context of SAIS models how and why $\alpha (i), \alpha _a(i), \delta _a(i)$ might depend on the actual prevalence i. While at low levels of prevalence the information received by susceptible hosts is unlikely to be confirmed by first-hand knowledge of actual cases, at moderate or high prevalence such confirmation by direct observation is more likely. Thus if such confirmation determines whether receiving information moves a given host into the A-compartment or the U-compartment, the ratio of newly created aware to newly created unwilling hosts may increase with the prevalence of the disease.

It is also of interest to investigate to what extent differences in the propensity to share information play a role in either inducing or enhancing sustained oscillations when the underlying disease dynamics is, for example, of type SIRS rather than SIS, or if demographics are included in the model.

Another possible direction of future research is to recast our models in the framework of stochastic processes. This might allow us to address the question of how long it takes to reach the absorbing disease-free state when populations are relatively small and trajectories are predicted to visit states with very low disease prevalence. It also would allow us to study the spread of the disease and of awareness on the relevant contact networks, similar to the work on nonreactive SAIS models in Sahneh et al. (2014).

Predictions of oscillatory dynamics have been reported for a number of previously studied models of behavioral epidemiology, both for models that incorporate demographics, and for models that, similarly to the ones studied here, ignore demographics. However, to the best of our knowledge, variability in the propensity of hosts to further disseminate awareness had not yet been identified as a possible driving mechanism for this phenomenon. We conclude this paper with a brief discussion of some related previous results. This review is not intended to be exhaustive; our goal is only to illustrate the variety of other mechanisms that appear to be capable of generating oscillations in disease prevalence.

Damped oscillations can be observed in non-reactive SAIS models when an asymptotically stable interior equilibrium has complex eigenvalues (Juher et al. 2015). They are also possible when the underlying disease dynamics is of type SIR without demographics. For example Epstein et al. (2008), considers spatially structured models where flight from endemic regions is identified as a key factor that will drive the overall dynamics.

Hopf bifurcations and sustained oscillations have been found in SEI- and SIR-based models with demographics (Lu et al. 2017) or with the related mechanism of prevalence-dependent recruitment of susceptibles into a core group (Velasco-Hernández et al. 1996). It is interesting to note that while the model of Lu et al. (2017) incorporates a reaction to media coverage, a type of awareness, it admits Hopf bifurcations even when the parameter that represents the strength of this reaction is zero. Sustained oscillations can also occur in models with demographics where the behavioral response represents compliance with vaccination and the driving mechanism involves real or perceived benefits of failing to adopt the behavioral response (Bauch 2005; d’Onofrio et al. 2007a, b; Reluga et al. 2006). Such benefits are absent in SAUIS models. They were also reported for models of vaccination compliance for influenza with fixed populations (Breban et al. 2007; Vardavas et al. 2007). In these models, vaccination decisions are made from year to year for different strains, so that the underlying disease dynamics is closer to type SIS than SIR. There is no direct information transfer in these models, and decisions are based on experience with past outbreaks rather than on current incidence levels. The latter features are also present in the models of d’Onofrio et al. (2007a, b).

Sustained oscillations also have been reported for a model that in some aspects resembles a reactive SAIS model, but makes an assumption about a fixed number of susceptibles. The authors of that paper conjectured that this assumption is needed for sustained oscillations in their model (see p. 1353 of Pawelek et al. 2014).

In Grassly et al. (2005) it is argued that oscillations observed in longitudinal studies of incidence data of syphilis may be entirely explainable by the underlying disease dynamics of type SIRS even though clear patterns of behavioral changes over time were present. This explanation is not consistent with an underlying disease dynamics of type SIS, and the authors of Grassly et al. (2005) partly base their conclusion on the absence of such oscillations in corresponding data sets on gonorrhea for the time frame of 1941–2001. However, as we mentioned in the introduction, a recent increase in the incidence of gonorrhea and other STIs that appears to be driven by changing behavior patterns has been reported (Wilton 2015).

Finally, there exists a substantial literature on oscillations in models with a behavioral response that involve rewiring the contact network (see, for example, Gross et al. 2006; Risau-Gusmán and Zanette 2009; Shaw and Schwartz 2008; Zhou et al. 2012). While these are individual-based models, analytical confirmation can sometimes be obtained by coarse-graining (Gross and Kevrekidis 2008) or using pair-approximations to build corresponding ODE-models with nonuniform mixing (Szabó-Solticzky et al. 2015; Szabó et al. 2012). Thus the mechanism that drives oscillations in these models is very different from the one in our SAUIS models, where uniform mixing is implicitly assumed.

Appendix: Proof of Lemma 2

Proof of Lemma 2

We can rewrite the Eq. (2) of the i-nullcline as

$$\begin{aligned} a_1(i) = \frac{\beta - \delta }{\beta - \beta _a} -\frac{\beta }{\beta - \beta _a}i. \end{aligned}$$

(17)

We will be using subscripts 1 and 2 in this proof to distinguish the function defined by (17) from the right-hand side of (4) that defines the positive branch $a_2(i)$ of the a-nullcline:

$$\begin{aligned} a_2(i)= & {} \frac{1}{2} \, \left( 1 - i - \frac{(\alpha _i(i) + \beta _a) i + \delta _a(i)}{\alpha _a(i)} \right. \nonumber \\&\left. + \ \sqrt{\left( 1 - i - \frac{(\alpha _i(i) + \beta _a) i + \delta _a(i)}{\alpha _a(i)} \right) ^2 + 4 \frac{\alpha _i(i)}{\alpha _a(i)} \, i \,(1-i) + \frac{4p(i) \delta }{\alpha _a(i)}i } \,\,\right) .\nonumber \\ \end{aligned}$$

(18)

Note that for $i = 1 - \frac{\delta }{\beta }$, we get

$$\begin{aligned} a_1\left( 1 - \frac{\delta }{\beta }\right) = 0 < a_2\left( 1 - \frac{\delta }{\beta }\right) . \end{aligned}$$

(19)

All biologically feasible interior equilibria must have their coordinate $i^*$ in the interval $\left( 0, 1 - \frac{\delta }{\beta }\right)$, as the graph of the linear function $a_1(i)$ cannot intersect $int(\Omega )$ for other values of i.

Note also that in this terminology Eq. (5) is equivalent to

$$\begin{aligned} a_1(0) > a_2 (0). \end{aligned}$$

(20)

Under the assumptions of part (a), this guarantees the existence of at least one $i^* \in \left( 0, 1 - \frac{\delta }{\beta }\right)$ with $a_1(i^*) = a_2(i^*)$. Moreover, the intersection of the graphs of $a_1(i), a_2(i)$ occurs in $int(\Omega )$. For $\delta \ge \beta _a$ this immediately follows because $a_1(i)$ lies entirely inside $\Omega$ for $i \in (0,1-\delta /\beta )$ since $a_1(0) < 1$.

For $\delta < \beta _a$, notice that $a_2(1) \ge 0$ and consider the values $i_0, i_1$ such that $a_2(i)$ enters $\Omega$ crossing the boundary $a+i=1$ at $(i_0,1-i_0)$, and $(i_1,1-i_1)$ is the intersection point of $a_1(i)$ with the same boundary. After substituting $a(i_0) = 1-i_0$ into the left-hand side of (3), solving for $i_0$, and simplifying we obtain:

$$\begin{aligned} i_0=\frac{\beta _a - \delta _a -p \delta + \sqrt{(\beta _a - \delta _a -p \delta )^2+4\beta _a\delta _a}}{2\beta _a} > i_1=1-\frac{\delta }{\beta _a}. \end{aligned}$$

Therefore, the graphs of $a_1(i)$ and $a_2(i)$ must intersect on the part of the graph of $a_1(i)$ that connects $\left( 1 - \frac{\delta }{\beta }, 0\right)$ and $(i_1,1-i_1)$, inside of $\Omega$.

To show uniqueness, under the assumptions of (a1), consider the function obtained by substituting the above expression for $a_1(i)$ into the left-hand side of (3). This substitution results in a quadratic function Q(i) that can have at most two roots. The combination of (19) and (20) immediately rules out the existence of $0< i^*< i^{**} < 1 - \frac{\delta }{\beta }$ at which the graphs of the functions $a_1(i)$ and $a_2(i)$ would actually cross. Moreover, if the the graphs were tangential at $i^\circ \in \{i^*, i^{**}\}$, then $i^\circ$ would be a root of Q of multiplicity 2, and the other equilibrium would need to be a third root of Q, which is impossible.

For the proof of part (b1), consider the limiting case where $\alpha _i = p = 0$. While we have explicitly ruled it out in our model assumptions, it will be useful in our argument. In this case, after substituting the expression for $a_1(i)$ in (17) into the left hand side of (3), we obtain a quadratic function in i:

$$\begin{aligned} H(i) = \left( \frac{\beta -\delta }{\beta - \beta _a} - \frac{\beta }{\beta -\beta _a}i \right) ^2 - \left( 1 - i - \frac{\beta _a i + \delta _a}{\alpha _a} \right) \left( \frac{\beta -\delta }{\beta -\beta _a} - \frac{\beta }{\beta - \beta _a}i \right) . \end{aligned}$$

Then there exists an interior equilibrium if and only if $H(i) = 0$ for some $0< i < 1 - \frac{\delta }{\beta }$. By directly evaluating H(i) at $i = 0$ and at $i = 1 - \frac{\delta }{\beta }$, we get

$$\begin{aligned} H(0) \le 0\ \ \ \text{ and }\ \ \ H\left( 1-\frac{\delta }{\beta }\right) = 0, \end{aligned}$$

(21)

where the first part relies on the assumption that inequality in (5) is reversed.

By evaluating $H'(i)$ at $i = 1 - \frac{\delta }{\beta }$, we get

$$\begin{aligned} H'\left( 1-\frac{\delta }{\beta }\right) = \frac{\alpha _a\delta + \beta _a\delta - \beta \beta _a - \beta \delta _a}{\alpha _a(\beta -\beta _a)}. \end{aligned}$$

Under the assumptions of our model, the denominator is positive, and inequality (6) is equivalent to $H'\left( 1-\frac{\delta }{\beta }\right) \ge 0$. It then follows from (21) that $H(i) < 0$ for all $0< i < 1-\frac{\delta }{\beta }$. Hence, the graphs of $a_1(i)$ and $a_2(i)$ do not intersect in $int(\Omega )$. Since the inequality (20) is assumed to be reversed, we can infer

$$\begin{aligned} a_2(i) > a_1(i)\ \ \ \text{ for } \text{ all }\ \ 0< i < 1 - \frac{\delta }{\beta }. \end{aligned}$$

(22)

When $\alpha _i + p > 0$, compared to the limiting case above where $\alpha _i = p = 0$, the expression for $a_1(i)$ stays the same while $a_2(i)$ becomes larger for each $0< i < 1 - \frac{\delta }{\beta }$. That is, in $int(\Omega )$, when $\alpha _i + p > 0$, we still have (22).

Therefore, there is no interior equilibrium.

For the proof of (a2), suppose $(i^*, a^*)$ is an interior equilibrium, which must be a point on the line segment $a_1(i)$ for some $0< i^* < 1 - \frac{\delta }{\beta }$. On this line segment, a(i) is decreasing, and $s(i) := 1 - a(i) - i$ is increasing, as can be seen if we rewrite (17) in its equivalent form given in (2) $i(a) = 1 - \frac{\delta }{\beta } - \left( 1-\frac{\beta _a}{\beta }\right) a$:

$$\begin{aligned} s(i) = 1 - a(i) - \left[ 1 - \frac{\delta }{\beta } - \left( 1-\frac{\beta _a}{\beta }\right) a(i) \right] = \frac{\delta }{\beta } - \frac{\beta _a}{\beta }a(i). \end{aligned}$$

(23)

We dropped the subscript here since we are only considering $a(i) = a_1(i)$ in this proof. Now we can rewrite the right hand side of the first line of (1), treating a and s as functions of i as in (17) and (23), to define the following function:

$$\begin{aligned} F(i)= & {} \alpha _i(i)s(i)i + \alpha _a(i)s(i)a(i) + p(i)\delta i - \beta _a a(i) i - \delta _a(i)a(i)\\= & {} \left( \alpha _i(i)\frac{s(i)i}{a(i)} + \alpha _a(i)s(i) - \beta _a i - \delta _a(i) + \frac{p(i)\delta i}{a(i)}\right) a(i)\\= & {} \left[ \alpha _i(i)\frac{\frac{\delta }{\beta }-\frac{\beta _a}{\beta }a(i)}{a(i)}i + \alpha _a(i)\left( \frac{\delta }{\beta } - \frac{\beta _a}{\beta }a(i)\right) - \beta _a i - \delta _a(i) + \frac{p(i)\delta i}{a(i)}\right] a(i). \end{aligned}$$

Let f(i) be the first factor so that $F(i) = f(i)a(i)$. It can be expressed as:

$$\begin{aligned} f(i)= & {} \alpha _i(i)g(i)i + \alpha _a(i)h(i) - \delta _a(i) + \frac{p(i)\delta i}{a(i)},\ \text{ where }\nonumber \\ g(i)= & {} \frac{\frac{\delta }{\beta } - \frac{\beta _a}{\beta }a(i)}{a(i)} = \frac{\delta }{\beta a(i)} - \frac{\beta _a}{\beta },\nonumber \\ h(i)= & {} \frac{\delta }{\beta } - \frac{\beta _a}{\beta }a(i) - \frac{\beta _a}{\alpha _a(i)}i\nonumber \\= & {} \frac{\delta }{\beta } - \beta _a\left( \frac{a(i)}{\beta } + \frac{i}{\alpha _a(i)}\right) \nonumber \\= & {} \frac{\delta }{\beta } - \beta _a\left( \frac{\frac{\beta -\delta }{\beta -\beta _a}-\frac{\beta }{\beta -\beta _a}i}{\beta } + \frac{i}{\alpha _a(i)}\right) \nonumber \\= & {} \frac{\delta }{\beta } - \frac{\beta _a}{\beta }\left( \frac{\beta -\delta }{\beta -\beta _a}-\frac{\beta }{\beta -\beta _a}i + \frac{\beta }{\alpha _a(i)}i\right) \nonumber \\= & {} \frac{\delta - \beta _a}{\beta - \beta _a} + \beta _a i\left( \frac{1}{\beta -\beta _a} - \frac{1}{\alpha _a(i)}\right) . \end{aligned}$$

(24)

Then at $(i^*, a^*)$ we must have $f(i^*) = 0$. Under the assumptions of (a2), g(i) and h(i) are increasing functions in i that take positive values at any interior equilibrium. Hence f(i) is a strictly increasing function in i and we can rule out the existence of two interior equilibria.

The result in part (b2) follows now by the same argument that was used in the proof of part (a1): Under the assumptions of (b2),

$$\begin{aligned} a_1\left( 1-\frac{\delta }{\beta }\right)< & {} a_2\left( 1 - \frac{\delta }{\beta }\right) ,\nonumber \\ a_1(0)< & {} a_2(0). \end{aligned}$$

Then an interior equilibrium could exist only if the graphs of $a_1(i)$ and $a_2(i)$ would cross at two or more of them, or if $a_1(i)$ and $a_2(i)$ are tangential at one such equilibrium. The former is ruled out by our uniqueness argument for part (a2). Note that this argument does not rely on inequality (5), which is only needed to prove existence. Situations where the graphs would cross at one equilibrium $i^*$ and be tangential at another equilibrium $i^{**}$ can also be ruled out on purely geometrical grounds: Note that $a_2(i)$ does not depend on $\beta$, but $a_1(i)$ does. When we slightly increase $\beta$, both $a_1(0)$ and $1 - \frac{\delta }{\beta }$ increase. Thus by slightly altering $\beta$ without violating the assumptions of (a2), we would obtain a parameter setting with three interior equilibria, which we have already shown to be impossible.

For the proof of (a3), consider again the function f(i) defined in (24) and rearrange its terms in the following way:

$$\begin{aligned} f(i)= & {} \alpha _i(i)\frac{\frac{\delta }{\beta } - \frac{\beta _a}{\beta }a(i)}{a(i)}i + \alpha _a(i)\left( \frac{\delta }{\beta } - \frac{\beta _a}{\beta }a(i)\right) - \beta _a i - \delta _a(i) + \frac{p(i)\delta i}{a(i)}\\= & {} \frac{\delta - \beta _a a(i)}{\beta }\left( \frac{\alpha _i(i)}{a(i)}i + \alpha _a(i)\right) + \frac{p(i)\delta i}{a(i)} - \beta _a i - \delta _a(i). \end{aligned}$$

Note that under the assumptions of (a3), the sum of the terms that enter this expression with a positive sign is strictly increasing in i, and the sum of the terms $-\beta _a i - \delta _a(i)$ is nondecreasing. Thus, f(i) is strictly increasing in i and we can rule out the existence of two interior equilibria.

The derivation of part (b3) from part (a3) is exactly analogous to the derivation of part (b2) from part (a2). $\square$

Remark

Note that the argument in the proof of (b1) relies only on the sign of $H'\left( 1-\frac{\delta }{\beta }\right)$. Additional sufficient conditions that preclude the existence of interior equilibria could be derived by considering other properties of H(i) or of its counterpart G(i) that allows for $\alpha _i + p > 0$. Our lemma is not meant to be exhaustive in this respect. For example, the leading coefficient of the quadratic function H(i) can be written as

$$\begin{aligned} \frac{\beta }{\beta - \beta _a}\left( \frac{\beta }{\beta - \beta _a} - 1 - \frac{\beta _a}{\alpha _a}\right) , \end{aligned}$$

and is positive if, and only if, $\alpha _a > \beta - \beta _a$, which together with (21) gives a proof of (b2) for the special case of a constant function $\alpha _a(i)$.