remarkRemark hypothesisHypothesis claimClaim Behaviour response to novel and
established diseasesA. Kaur, R. C. Tyson, and I. R. Moyles
[][nocite]ex_supplement

# The impact of fear and behaviour response to established and novel diseases
††thanks: Submitted to the editors DATE. This work was funded by the Natural
Sciences and Engineering Research Council of Canada through the Discovery
Grant and Emerging and Infectious Disease Modelling programs.

Avneet Kaur Department of Computer Science, Mathematics, Physics, and
Statistics, University of British Columbia, Kelowna, ON (, ).
<EMAIL_ADDRESS><EMAIL_ADDRESS>Rebecca C. Tyson
22footnotemark: 2 Iain R. Moyles Department of Mathematics and Statistics,
York University, Toronto, ON (, Corresponding Author<EMAIL_ADDRESS>
###### Abstract

We analyze a disease transmission model that allows individuals to acquire
fear and change their behaviour to reduce transmission. Fear is acquired
through contact with infected individuals and through the influence of fearful
individuals. We analyze the model in two limits: First, an Established Disease
Limit (EDL), where the spread of the disease is much faster than the spread of
fear, and second, a Novel Disease Limit (NDL), where the spread of the disease
is comparable to that of fear. For the EDL, we show that the relative rate of
fear acquisition to disease transmission controls the size of the fearful
population at the end of a disease outbreak, and that the fear-induced contact
reduction behaviour has very little impact on disease burden. Conversely, we
show that in the NDL, disease burden can be controlled by fear-induced
behaviour depending on the rate of fear loss. Specifically, fear-induced
behaviour introduces a contact parameter $p$, which if too large prevents the
contact reduction from effectively managing the epidemic. We analytically
identify a critical prophylactic behaviour parameter $p=p_{c}$ where this
happens leading to a discontinuity in epidemic prevalence. We show that this
change in disease burden introduces multiple epidemic waves.

###### keywords:

disease model; behaviour response; multiple outbreaks; disease bifurcation;
fear.

92D30;35B40;34C23;34E20

##  1 Introduction

Human behavioural changes play a crucial role in shaping the trajectory of
infectious diseases. As the disease pathogen spreads, the associated morbidity
and mortality give rise to population-level perceptions of the risk and
severity of the disease. As a result, people try to minimise those risks by
spontaneously changing their behaviours, even in the absence of government
directives [2, 3].

Most commonly, risk minimisation is done through prophylactic actions such as
social distancing, wearing face masks, etc. Modelling behaviour and opinion
dynamics alone can lead to rich patterns in behaviour adoption and complex
dynamics such as the degree to which opinion strength affects behaviour [4,
5]. However, when coupled with disease dynamics, behaviour can shape the
course of the epidemic. This observation has led to the important field of
Behavioural Epidemiology which emphasizes the incorporation of theories from
sociology and psychology as key components disease models [1, 19]. Several
reviews including behaviour and its importance in infectious disease modelling
have been written [16, 38].

Many modelling efforts in Behavioural Epidemiology have demonstrated that
spontaneous behavioural responses to increased disease prevalence can result
in disease dynamics that differ substantially from those predicted by models
that do not include human behaviour. For example, incorporating prophylactic
behaviour in the response to an emerging disease can induce multiple waves of
infection [28, 37, 22]. Furthermore, including behaviour in disease models has
produced mechanisms that could explain trends observed in data. For example, a
model of individual health decisions captured face mask usage during the 2003
SARS epidemic in Hong Kong [13], a model of risk perception was able to
explain a sudden change in disease incidence for the 2009 H1N1 epidemic in
Italy [35], and a model of perceived immunity was able to explain the fall
2021 wave of COVID-19 in Ontario Canada [27]. It has even been demonstrated
that incorporating the planning horizon time into behavioural decision making
can itself impact disease burden [31].

Naturally, as the inclusion of behaviour strengthened epidemiological
predictions, the structure of these mathematical models and the mechanisms
that affect disease burden became an area of interest. For example, an
adaptation of the classic Susceptible-Infected-Recovered (SIR) formulation of
Kermack and McKendrick [23] introduced fear as a second contagion and showed
that this dual-contagion model can lead to multiple epidemic waves [15]. Perra
et al (2011) generalized this model framework and showed that when fear
spreads at a rate faster than disease, the fear-induced contact reduction
reduces the final epidemic size (i.e., the total infected population) [32] . A
similar bifurcation in final epidemic size was observed in Poletti et al
(2011) [35] where they showed that behaviour changes driven by a high
perception of risk of disease, fast adoption of fear, and a relatively slow
fear loss, can significantly change the disease dynamics, particularly sudden
changes from slow to fast growth in incidence case numbers.

The complexity of the behaviourally-mediated disease dynamics has limited
model analysis. Perra et al (2011) generalized some fear-responsive models and
used linear theory to identify some bifurcations in final sizes of fear
compartments [32]. However, their most interesting observation was a
discontinuity in epidemic prevalence when a self-reinforcement mechanism of
fear loss was included in their model. The authors concluded that the
mechanism of this bifurcation was outside of linear theory and left it
unresolved. A similar fear-loss mechanism was introduced by Epstein et al
where they observed bifurcations in epidemic prevalence, but through
simulations only [14]. Yet another bifurcation in final epidemic size was
observed in the qualitative analysis of a slightly different model [33]. This
model included imitative fear loss, which corresponds to a fearful population
losing fear upon interaction with regular susceptibles ($S$), who are
unafraid. Overall, the relationship between different disease dynamics,
particularly the bifurcations that separate them, and the underlying
parameters that drive them is key ingredient in developing our understanding
of epidemics and for improving public health interventions.

We consider a fear-incorporating disease model under two disease scenarios,
one where the circulating disease is established (e.g., the seasonal flu) and
one where the circulating disease is novel (e.g., COVID-19 in early 2020).
Using asymptotic and weakly non-linear analysis we determine the role of
behaviour in shaping disease burden outcomes for each scenario, including the
resolution of discontinuities in epidemic prevalence first observed in [32].

The outline of our article is as follows. In Section 2 we define and non-
dimensionalize the fear-disease model, determine the established and novel
disease limits, and discuss equilibria and stability. We analyse the
Established Disease Limit in Section 3 and the Novel Disease Limit in Section
4. We discuss the implications of our results and offer conclusions in Section
5.

##  2 Model description

The model we consider has developed from a series of papers (cf. [15, 32,
14]); we consider the version from Epstein et al (2011) [14]. The model is an
extension of the fundamental Susceptible-Infected-Recovered model [23] in
which the susceptible compartment is split into two groups: the susceptibles
who are fearless ($S$) and fearful ($S_{\rm{fd}}$) of the disease. The
compartment $I$ refers to the infectious individuals, and $R_{\rm{nat}}$
refers to the individuals who recovered naturally from the disease. The
compartmental diagram of the model is shown in Fig. 1.

Figure 1: Compartmental diagram of the model [14]. The compartments $S$ and
$S_{\rm{fd}}$ refer to the fearless and disease-fearful susceptibles,
respectively. The compartment $I$ represents the infectious individuals, and
the compartment $R_{\rm{nat}}$ refers to the individuals who recovered
naturally from the disease. The various transition rates between compartments
appear above each arrow.

We assume that each state variable, $X\in[S,S_{\rm{fd}},I,R_{\rm{nat}}]$, in
our model represents the corresponding fraction of the total population, hence
$X\in[0,1]$ and are subject to the constant population constraint
$S+S_{\rm{fd}}+I+R_{\rm{nat}}=1$. Assuming a well-mixed population, the
differential equations corresponding to Fig. 1 are

(1a) $\displaystyle\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!S}{\mathop{\mathrm{\mathstrut d}}\nolimits\\!t}$
$\displaystyle=-\beta
IS-\beta_{\rm{fd}}(S_{\rm{fd}}+I)S+\gamma_{\rm{f}}S_{\rm{fd}}+\alpha_{\rm{f}}(1-S-S_{\rm{fd}}-I)S_{\rm{fd}},$
(1b) $\displaystyle\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!S_{\rm{fd}}}{\mathop{\mathrm{\mathstrut d}}\nolimits\\!t}$
$\displaystyle=-p\beta
IS_{\rm{fd}}+\beta_{\rm{fd}}(S_{\rm{fd}}+I)S-\gamma_{\rm{f}}S_{\rm{fd}}-\alpha_{\rm{f}}(1-S-S_{\rm{fd}}-I)S_{\rm{fd}},$
(1c) $\displaystyle\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!I}{\mathop{\mathrm{\mathstrut d}}\nolimits\\!t}$
$\displaystyle=\beta IS+p\beta IS_{\rm{fd}}-\gamma I,$ (1d)
$\displaystyle\frac{\mathop{\mathrm{\mathstrut
d}}\nolimits\\!R_{\rm{nat}}}{\mathop{\mathrm{\mathstrut d}}\nolimits\\!t}$
$\displaystyle=\gamma I.$

where $\beta$ is the disease transmission rate, $\beta_{\rm fd}$ the fear
transmission rate, $\gamma$ the disease recovery rate, $\gamma_{\rm f}$ the
rate of spontaneous loss of fear, $\alpha_{\rm f}$ the rate of loss of fear
through contact, and $p$ the relative reduction in rate of disease
transmission

The role of fear impacts the model Eq. 1 in four ways. There are two
mechanisms of fear acquisition which come from susceptible individuals
interacting with either infected individuals or those susceptible but fearful.
There are also two mechanisms of fear loss with a natural loss through which
fearful susceptibles lose their fear and a self-reinforcing loss from fearful
susceptible individuals who interact with those who have recovered from the
disease. We consider this second mechanism to be fear loss due to complacency,
as the observation of those who have recovered from the disease leads to a
lower perception of risk and reduced prophylactic behaviour. It is tuning of
this mechanism that yielded discontinuous epidemic prevalence in other studies
[32, 14]. Numerical study of the model Eq. 1 demonstrated that the final
epidemic size increases suddenly as either $p$ or $\alpha_{\rm{f}}$ increase
past a critical value [14]. This same study also showed that, as $p$
increases, i.e., as fear-induced contact reduction becomes less effective, the
infection dynamics transition from a pair of waves to a single small infection
wave.

###  2.1 Steady State Analysis

We review some of the rich aspects of the phase space of Eq. 1 as first noted
by Perra et al (2011) [32]. Particularly, we note that Eq. 1 permits two
classes of equilibria,

(2a) $\displaystyle E_{0}(S^{*})=$ $\displaystyle(S^{*},0,0,1-S^{*}),$ (2b)
$\displaystyle E_{1}(S_{\rm{fd}}^{*})=$
$\displaystyle(\bar{S},S_{\rm{fd}}^{*},0,1-\bar{S}-S_{\rm{fd}}^{*}),$

where,

(3)
$\bar{S}=\frac{\gamma_{\rm{f}}+\alpha_{\rm{f}}(1-S_{\rm{fd}}^{*})}{\alpha_{\rm{f}}+\beta_{\rm{fd}}},$

and $S^{*}$ and $S_{\rm{fd}}^{*}$ are the fractions of the population in the
regular susceptible group ($S$), and the fearful susceptible group
($S_{\rm{fd}}$) at equilibrium. As is often the case in disease models, these
values cannot be determined uniquely via steady-state analysis as they are
dependent on the initial conditions of the model.

The equilibrium $E_{0}$ Eq. 2a has no fearful susceptible individuals, no
infected individuals, and a recovered population of $1-S_{0}^{*}$. When
$S_{0}^{*}=1$ then there is no recovered population and the equilibrium
$E_{0}(1)$ is classically known as a disease-free equilibrium (DFE) [11] since
the disease never spread. When $S_{0}^{*}\neq 0$ then this equilibrium
represents a population that has recovered from the disease and retains no
fear thereof. We therefore call equilibrium Eq. 2a the fear-free equilibrium.

The second equilibrium, $E_{1}$ (2b), differs from $E_{0}$ in that it has a
non-zero fearful susceptible population. We therefore call this equilibrium
the fear-endemic equilibrium. We note that it is possible to have a fear-
endemic disease-free equilibrium meaning that fear spreads even when no
disease is present. This occurs if $S_{\rm{fd}}^{*}+\bar{S}=1$ as there are no
recovered individuals in the system indicating that the disease did not
spread.

####  2.1.1 Stability of the Disease-Free Equilibrium

To investigate the stability of the DFE $E_{0}(1)=(1,0,0,0)$, we obtain the
eigenvalues of the Jacobian around the DFE. The Jacobian has a zero eigenvalue
and two non-zero eigenvalues. The zero eigenvalue is related to the degenerate
$E_{0}(S^{*})$ for multiple values of $S^{*}$ and does not affect stability of
the DFE to large perturbations. Hence, we consider the non-zero eigenvalues of
the Jacobian

(4a) $\displaystyle\lambda_{1}$ $\displaystyle=\beta-\gamma,$ (4b)
$\displaystyle\lambda_{2}$ $\displaystyle=\beta_{\rm{fd}}-\gamma_{\rm{f}}.$

From here we can define a disease basic reproduction number $\mathcal{R}_{0}$
and a fear basic reproduction number $\mathcal{R}_{\rm f}$,

(5) $\displaystyle\mathcal{R}_{0}=\frac{\beta}{\gamma},\quad\mathcal{R}_{\rm
f}=\frac{\beta_{\rm fd}}{\gamma_{\rm f}},$

such that the DFE is unstable if either $\mathcal{R}_{0}>1$ or
$\mathcal{R}_{\rm f}>1$. These reproduction numbers could also be derived
using the next generation method [12].

We are interested in understanding the dynamics of the system when both the
disease and fear contagions can spread in the population, so we assume that
both reproduction numbers in Eq. 5 are greater than one so that the
eigenvalues Eq. 4 are positive.

####  2.1.2 Stability of the Fear-Free Equilibrium

To determine the stability of the fear-free equilibrium when the disease has
spread, i.e., $E_{0}(S^{*})$, $S^{*}<1$ in Eq. 2a, we obtain the non-zero
eigenvalues of the Jacobian at $E_{0}(S^{*})$:

(6a) $\displaystyle\lambda_{1}$
$\displaystyle=(S^{*}\mathcal{R}_{0}-1)\gamma,$ (6b)
$\displaystyle\lambda_{2}$ $\displaystyle=(S^{*}\mathcal{R}_{\rm
f}-1)\gamma_{\rm{f}}-(1-S^{*})\alpha_{\rm{f}}.$

Both of these eigenvalues are negative if

(7) $\displaystyle
S^{*}<\min\left(\frac{1}{\mathcal{R}_{0}},\frac{a+1}{a+\mathcal{R}_{\rm
f}}\right);\qquad a=\frac{\alpha_{\rm{f}}}{\gamma_{\rm{f}}}.$

In the absence of fear loss due to complacency (when $a=0$), condition Eq. 7
reduces to

$\max(\mathcal{R}_{0}^{\rm eff},\mathcal{R}_{\rm f}^{\rm eff})<1,$

where $\mathcal{R}_{i}^{\rm eff}=S^{*}\mathcal{R}_{i}$ are the effective
reproduction numbers of the final fearless susceptible population $S^{*}$ and
appear regularly in disease models as a condition of reaching herd immunity
[17, 18].

####  2.1.3 Stability of the Fear-Endemic Equilibrium

The non-zero eigenvalues of the Jacobian at the fear-endemic equilibrium,
$E_{1}(S_{\rm{fd}}^{*})$ from Eq. 2b, are:

(8a) $\displaystyle\lambda_{1}$ $\displaystyle=-\mathcal{R}_{\rm
f}\gamma_{\rm{f}}S_{\rm{fd}}^{*},$ (8b) $\displaystyle\lambda_{2}$
$\displaystyle=\left[\mathcal{R}_{0}S_{\rm{fd}}^{*}\left(p-\frac{a}{a+\mathcal{R}_{\rm
f}}\right)+\mathcal{R}_{0}\left(\frac{a+1}{a+\mathcal{R}_{\rm
f}}\right)-1\right]\gamma.$

The first of these eigenvalues is always negative, but there are a myriad of
stability conditions for Eq. 8b. In the absence of fear loss due to
complacency ($a=0$) then Eq. 8b becomes

(9)
$\displaystyle\lambda_{2}=\left[\mathcal{R}_{0}S_{\rm{fd}}^{*}p+\frac{\mathcal{R}_{0}}{\mathcal{R}_{\rm
f}}-1\right]\gamma.$

The first term, $\mathcal{R}_{0}S_{\rm{fd}}^{*}p$ is another effective
reproduction number where $p$, the modified contact parameter, reduces the
disease transmissibilty for susceptible people who are fearful. From Eq. 9 we
immediately see that if $\mathcal{R}_{0}>\mathcal{R}_{\rm f}$ then the fear-
endemic equilibrium $E_{1}(S_{\rm{fd}}^{*})$ is unstable. This result is
intuitive because, in this case, the growth rate of the disease is larger than
that of fear. However, even if $\mathcal{R}_{\rm f}>\mathcal{R}_{0}$, the
fear-endemic equilibrium can still be unstable if the behaviour-modified
effective reproduction number $\mathcal{R}_{0}S_{\rm{fd}}^{*}p$ is large
enough, which happens if the fearful group has only a weak prophylactic
response (large $p$) and/or a large population size. Thus, even if fear grows
quickly, if the prophylactic behaviour that is generated from that fear is not
strong enough, the equilibrium can become unstable.

Overall, the stability criteria for the steady states clearly demonstrates
that rich dynamics are possible and the non-linear responses observed by Perra
et al (2011) are unsurprising [32]. This motivates the deeper non-linear
analysis of these models which we now perform.

###  2.2 Model Scaling

The steady state analysis above uncovered key parameter groupings,
particularly the reproduction numbers $\mathcal{R}_{i}$ Eq. 5 and the fear-
loss ratio, $a$ Eq. 7. We are therefore motivated to non-dimensionalize the
model Eq. 1. We scale time $t\sim\gamma^{-1}$ with the disease recovery time
to produce the non-dimensional model

(10a) $\displaystyle\dot{S}=$
$\displaystyle-\mathcal{R}_{0}SI-\delta\left[\mathcal{R}_{\rm
f}(S_{\rm{fd}}+I)S-S_{\rm{fd}}-a(1-S-S_{\rm{fd}}-I)S_{\rm{fd}}\right],$ (10b)
$\displaystyle\dot{S}_{\rm fd}=$
$\displaystyle-\mathcal{R}_{0}pS_{\rm{fd}}I+\delta\left[\mathcal{R}_{\rm
f}(S_{\rm{fd}}+I)S-S_{\rm{fd}}-a(1-S-S_{\rm{fd}}-I)S_{\rm{fd}}\right],$ (10c)
$\displaystyle\dot{I}=$
$\displaystyle\mathcal{R}_{0}\left(S+pS_{\rm{fd}}\right)I-I,$

where the over-dot indicates differentiation with respect to non-dimensional
time. The reproduction numbers and fear-loss ratio appear in this model along
with a new non-dimensional parameter

(11) $\displaystyle\delta=\frac{\gamma_{\rm f}}{\gamma},$

which is the ratio of the rates of recovery for fear and disease. We make the
assumption that $\delta\ll 1$ meaning that the average natural recovery time
for fear is much longer than that of the disease. This assumption is supported
through literature where studies on fear extinction have indicated that
therapy is often required to eliminate fear and that there can be high rates
of relapse [9, 36, 30]. Conversely, a study of elderly respondents with a
common cold showed that only 15% had not recovered within 8 days [24].
Similarly, a study on COVID-19 found that 80% of people had recovered within
one month [26]. The recovery rate $\gamma$ in our model represents the
timescale over which someone is no longer infectious which is likely faster
than the time it takes for a patient to self-declare as being recovered.
Recovery studies thus likely overestimate the loss of infectiousness time
(recovery time in our context).

For the non-dimensional model we take as our initial condition a susceptible
population that is entirely fearless toward the disease, and a small group of
infected individuals. All other subgroups have zero population. Thus, we take

(12) $\displaystyle S(0)=1-\delta I_{0},\quad S_{\rm{fd}}(0)=0,\quad
I(0)=\delta I_{0},$

noting that there is no loss of generality if we allow the initial population
to include some disease-fearful individuals, i.e. that $S_{\rm{fd}}(0)>0$.
Based on the initial conditions Eq. 12, it is natural to transform the
variables as

(13) $\displaystyle S=1+\delta x,\quad S_{\rm{fd}}=\delta y,\quad I=\delta z,$

leading to the new problem

(14a) $\displaystyle\dot{x}=$ $\displaystyle-\mathcal{R}_{0}z(1+\delta
x)-\delta\left[\mathcal{R}_{\rm f}(y+z)(1+\delta x)-y-a(x+y+z)y\right],$ (14b)
$\displaystyle\dot{y}=$ $\displaystyle-\delta\mathcal{R}_{\rm
f}pyz+\delta\left[\mathcal{R}_{\rm f}(y+z)(1+\delta x)-y-a(x+y+z)y\right],$
(14c) $\displaystyle\dot{z}=$ $\displaystyle\mathcal{R}_{0}(1+\delta x+\delta
py)z-z.$

When all non-dimensional parameters are $\mathcal{O}(1)$ then since $\delta\ll
1$ the leading order problem is

(15a) $\displaystyle\dot{x}=$ $\displaystyle-\mathcal{R}_{0}z,$ (15b)
$\displaystyle\dot{y}=$ $\displaystyle 0,$ (15c) $\displaystyle\dot{z}=$
$\displaystyle(\mathcal{R}_{0}-1)z,$

and so we see that there is minimal growth in the fearful susceptible
population and that susceptibles quickly become infected,

(16) $\displaystyle z(t)=I_{0}\textrm{e}^{{(\mathcal{R}_{0}-1)t}}$

with an outbreak occurring under the assumption that $\mathcal{R}_{0}>1$. In
this limit the fear-reproduction number is comparable to the disease
reproduction number meaning that fear acquisition does not occur exceptionally
quickly. As such we refer to this as the Established Disease Limit (EDL).

Fear will grow slowly in Eq. 14 unless $\mathcal{R}_{\rm f}\gg 1$ indicating
that there is a rapid acquisition of fear. In this case we can write,
$\mathcal{R}_{\rm f}=\delta^{-1}\hat{\mathcal{R}}_{\rm f}$ with
$\hat{\mathcal{R}}_{\rm f}=\beta_{\rm fd}\gamma^{-1}$. Similarly scaling,
$a=\delta^{-1}\hat{a}$ with $\hat{a}=\alpha_{\rm f}\gamma^{-1}$ changes the
leading order of Eq. 14 to

(17a) $\displaystyle\dot{x}=$
$\displaystyle-\mathcal{R}_{0}z-\hat{\mathcal{R}}_{\rm f}(y+z),$ (17b)
$\displaystyle\dot{y}=$ $\displaystyle\hat{\mathcal{R}}_{\rm f}(y+z),$ (17c)
$\displaystyle\dot{z}=$ $\displaystyle(\mathcal{R}_{0}-1)z.$

Once again the disease grows exponentially via Eq. 16, but now the fearful
susceptible population grows with solution

(18) $\displaystyle y(t)=\frac{I_{0}\hat{\mathcal{R}}_{\rm
f}}{1-\mathcal{R}_{0}+\hat{\mathcal{R}}_{\rm
f}}\left(\textrm{e}^{{\hat{\mathcal{R}}_{\rm
f}t}}-\textrm{e}^{{(\mathcal{R}_{0}-1)t}}\right),$

and people become rapidly afraid. In fact, if $\hat{\mathcal{R}}_{\rm
f}>\mathcal{R}_{0}-1$ then they become afraid more rapidly than the disease
grows. We call the scenario where $\mathcal{R}_{\rm
f}\sim\mathcal{O}(\delta^{-1})$ the Novel Disease Limit (NDL) since fear
acquisition is very rapid suggesting that there is limited population exposure
to the pathogen. An example of such a disease would be COVID-19 in early 2020
before it approached endemicity. We analyze these two limits separately and
investigate the impact of fear and behaviour control on disease outcomes.

##  3 The Established Disease Limit (EDL)

As was outlined in Section 2.2 the EDL is one in which each reproduction
number leads to an instability in the disease and fear-free equilibrium
$\mathcal{R}_{i}>1$, but that $\mathcal{R}_{i}\sim\mathcal{O}(1)$. In this
case the fear grows slowly while the disease grows rapidly and therefore to
analyze the model we consider the initial conditions $S(0)=1-I_{0}$,
$S_{\rm{fd}}(0)=0$, and $I(0)=I_{0}$ noting that allowing an initial fearful
susceptible population does not alter the fundamental structure of the results
that follow. We take $S=x$, $S_{\rm{fd}}=\delta y$, and $I=z$ and write Eq. 10
as

(19a) $\displaystyle\dot{x}=$
$\displaystyle-\mathcal{R}_{0}xz-\delta\left[\mathcal{R}_{\rm f}(\delta
y+z)x-\delta y-\delta a(1-x-\delta y-z)y\right],$ (19b)
$\displaystyle\dot{y}=$ $\displaystyle-\mathcal{R}_{0}pyz+\mathcal{R}_{\rm
f}(\delta y+z)x-\delta y-\delta a(1-x-\delta y-z)y,$ (19c)
$\displaystyle\dot{z}=$ $\displaystyle\mathcal{R}_{0}(x+\delta py)z-z.$

We form an asymptotic expansion,

(20) $\displaystyle x\sim x_{0}+\delta x_{1}+\dots,\quad y\sim y_{0}+\delta
y_{1}+\dots,\quad z\sim z_{0}+\delta z_{1}+\dots,$

and to leading order Eq. 19 becomes

(21a) $\displaystyle\dot{x}_{0}=$ $\displaystyle-\mathcal{R}_{0}x_{0}z_{0}$
(21b) $\displaystyle\dot{y}_{0}=$ $\displaystyle\mathcal{R}_{\rm
f}x_{0}z_{0}-\mathcal{R}_{0}py_{0}z_{0},$ (21c) $\displaystyle\dot{z}_{0}=$
$\displaystyle(\mathcal{R}_{0}x_{0}-1)z_{0},$

where we see that the fearful susceptible population has decoupled from the
susceptible and infected group. Because of this, the leading order disease
dynamics follow a traditional SIR model which has an analytic solution (see
for example [7]),

(22) $\displaystyle
z(x)=1-x_{0}+\frac{1}{\mathcal{R}_{0}}\log\left(\frac{x_{0}}{1-I_{0}}\right),$

accompanied by the classic SIR results such as single-peak epidemic dynamics.
The leading order fearful susceptible population is determined by solving Eq.
21b yielding,

(23) $\displaystyle y_{0}=\frac{\mathcal{R}_{\rm
f}}{\mathcal{R}_{0}(1-p)}x_{0}^{p}\left((1-I_{0})^{1-p}-x_{0}^{1-p}\right).$

We note that Eq. 22 only provides an implicit relationship between $x_{0}$ and
$z_{0}$, but this can be used to furnish integral solutions explicit on time.
However, since the dynamics mimic the SIR model, the main outstanding issue is
the final size of the fearful susceptible group since we showed in Section 2
that the equilibria permit both fearless and fearful possibilities. As such,
we are primarily concerned with the final size of each compartment. There is
no endemic state from the disease due to an absence of waning immunity and so
from both Eq. 22 and the equilibrium discussion in Section 2 then
$z_{0}^{\infty}=\lim_{t\to\infty}z_{0}=0$. We can add Eq. 21a and Eq. 21c to
determine that the final size for the susceptible group.
$x_{0}^{\infty}=\lim_{t\to\infty}x_{0}$, satisfies,

(24)
$\displaystyle\frac{1}{\mathcal{R}_{0}}\log\left(\frac{1-I_{0}}{x_{0}^{\infty}}\right)=1-x_{0}^{\infty},$

and then we can use this value in Eq. 23 to determine $y_{0}^{\infty}$. We
thus far have that at steady state $S\sim x_{0}^{\infty}$ while
$S_{\rm{fd}}\sim\delta y_{0}^{\infty}$ and we therefore require
$x_{1}^{\infty}$ so that both final sizes are known to the same order. The
$\mathcal{O}(\delta)$ problem from Eq. 19 is

(25a) $\displaystyle\dot{x}_{1}=$
$\displaystyle-\mathcal{R}_{0}(x_{1}z_{0}+x_{0}z_{1})-\mathcal{R}_{\rm
f}x_{0}z_{0},$ (25b) $\displaystyle\dot{y}_{1}=$
$\displaystyle-\mathcal{R}_{0}p(y_{1}z_{0}+y_{0}z_{1})+\mathcal{R}_{\rm
f}(x_{1}z_{0}+x_{0}z_{1}+x_{0}y_{0})-y_{0}-a(1-x_{0}-z_{0})y_{0},$ (25c)
$\displaystyle\dot{z}_{1}=$
$\displaystyle(\mathcal{R}_{0}x_{0}-1)z_{1}+\mathcal{R}_{0}x_{1}z_{0}+\mathcal{R}_{0}py_{0}z_{0},$

subject to $x_{1}(0)=y_{1}(0)=z_{1}(0)=0$. Adding Eqs. 25a and 25c and
integrating to infinite time yields,

(26) $\displaystyle
x_{1}^{\infty}+y_{0}^{\infty}=-\int_{0}^{\infty}z_{1}\mathop{\mathrm{\mathstrut
d}}\nolimits\\!t.$

If we write $x_{1}=x_{0}u$ then we can transform Eq. 25a to

$\dot{u}=-\mathcal{R}_{0}z_{1}-\mathcal{R}_{\rm f}z_{0},\qquad u(0)=0,$

which upon integrating over infinite time yields

$u^{\infty}=\mathcal{R}_{0}(x_{1}^{\infty}+y_{0}^{\infty})+\mathcal{R}_{\rm
f}(x_{0}^{\infty}-1),$

where we have used that

$\int_{0}^{\infty}z_{0}=1-x_{0}^{\infty}.$

Finally then we have that

(27) $\displaystyle
x_{1}^{\infty}=\left(\frac{\mathcal{R}_{0}y_{0}^{\infty}-\mathcal{R}_{\rm
f}(1-x_{0}^{\infty})}{1-\mathcal{R}_{0}x_{0}^{\infty}}\right)x_{0}^{\infty}.$

###  3.1 Long-time limit

The final-size prediction for the fearful susceptible population is
$S_{\rm{fd}}\sim\delta y_{0}$ with $y_{0}$ given by Eq. 23. This is positive
which implies that there are always fearful susceptibles, something we know
not to be the case from the equilibrium considerations in Section 2. It turns
out (see Appendix A) that $S_{\rm{fd}}$ grows without bound and thus
eventually violates the assumption that $S_{\rm{fd}}\sim\mathcal{O}(\delta)$.
This happens in a timescale $t\sim\mathcal{O}(\delta^{-1})$ and so it suggests
that we let $t=\delta^{-1}\tau$ and write Eq. 1 in this new timescale as,

(28a) $\displaystyle S_{\tau}=$ $\displaystyle
S_{\rm{fd}}+a(1-S-S_{\rm{fd}})S_{\rm{fd}}-\mathcal{R}_{\rm f}S_{\rm{fd}}S,$
(28b) $\displaystyle{S_{\rm{fd}}}_{\tau}=$ $\displaystyle\mathcal{R}_{\rm
f}S_{\rm{fd}}S-S_{\rm{fd}}-a(1-S-S_{\rm{fd}})S_{\rm{fd}},$

with initial conditions $S(0)=x_{0}^{\infty}+\delta x_{1}^{\infty}$ and
$S_{\rm{fd}}(0)=\delta y_{0}^{\infty}$. In writing Eq. 28 we have used
extinction of the disease to eliminate $I$ in the long-time limit (see
Appendix A). Thus, fear remains the only contagious agent.

We add the component equations in Eq. 28, integrate, and use the initial
conditions to conclude that

(29) $\displaystyle
S+S_{\rm{fd}}=\eta;\qquad\eta=x_{0}^{\infty}+\delta(x_{1}^{\infty}+y_{0}^{\infty}).$

From this we can eliminate $S$ in Eq. 28b to determine the long-time fearful
susceptible population,

(30) $\displaystyle S_{\rm{fd}}(\tau)=\frac{(\kappa-1)\delta
y_{0}}{\delta\mathcal{R}_{\rm f}y_{0}-(\mathcal{R}_{\rm f}\delta
y_{0}+1-\kappa)\textrm{e}^{{(1-\kappa)\tau}}};\qquad\kappa=\mathcal{R}_{\rm
f}\eta-a(1-\eta),$

If $\kappa<1$ then the exponential in the denominator of Eq. 30 grows as
$\tau\to\infty$ so $S_{\rm{fd}}\to 0$. Conversely, if $\kappa>1$ then the
exponential decays and

(31) $\displaystyle S_{\rm{fd}}^{\infty}=\frac{\kappa-1}{\mathcal{R}_{\rm
f}}.$

From this we have that

(32) $\displaystyle S^{\infty}=\eta-
S_{\rm{fd}}^{\infty}=\frac{1+a(1-\eta)}{\mathcal{R}_{\rm f}}$

which aligns with the predicted steady state value Eq. 2b in Section 2 when
$S_{\rm{fd}}^{*}=S_{\rm{fd}}^{\infty}$. Thus the final fear size is dependent
on a new effective reproduction number $\kappa$. Since $\eta$ is the final
size of the two susceptible groups during the disease outbreak then $1-\eta$
is the population of those who have recovered. Therefore, fear can spread if
the fear reproduction number is sufficiently high such that infecting the
remaining susceptibles, $\eta$ with fear is greater than the fear extinction
from recovery $a(1-\eta)$. Since $\eta$ and $\mathcal{R}_{\rm f}$ are
independent of $a$ then a fear-endemic state, $S_{\rm{fd}}>0$ occurs when

(33) $\displaystyle a<a_{c};\qquad a_{c}=\frac{\mathcal{R}_{\rm
f}\eta-1}{1-\eta}.$

we demonstrate this in Fig. 2a for parameters $\mathcal{R}_{\rm f}=8$,
$\mathcal{R}_{0}=2$, $p=0.25$, and $I_{0}=\delta=0.01$ whence $\eta=0.2197$
and $a_{c}=0.9708$. If instead we fix $a$ then $\eta=\eta(\mathcal{R}_{\rm
f})$ and so the bifurcation relationship is non-linear. The bifurcation point,
${\mathcal{R}_{\rm f}}_{c}$ for a fear endemic equilibrium in this case from
Eq. 30 satisfies

(34) $\displaystyle{\mathcal{R}_{\rm f}}_{c}\eta({\mathcal{R}_{\rm
f}}_{c})-a(1-\eta({\mathcal{R}_{\rm f}}_{c}))=1$

and since $\delta\ll 1$ we pose an expansion ${\mathcal{R}_{\rm
f}}_{c}={\mathcal{R}_{\rm f}}_{c0}+\delta{\mathcal{R}_{\rm f}}_{c1}+\dots$
which yields,

(35) $\displaystyle{\mathcal{R}_{\rm
f}}_{c0}=\frac{1+a(1-x_{0})}{x_{0}},\quad{\mathcal{R}_{\rm
f}}_{c1}=-\frac{(1+a)(1+a(1-x_{0}))(y_{0}+x_{1})}{x_{0}^{3}}.$

For parameters $a=0$, $\mathcal{R}_{0}=2$, $p=0.25$, and $I_{0}=\delta=0.01$
we compute that ${\mathcal{R}_{\rm f}}_{c}=4.693$ which is demonstrated in
Fig. 2b to good agreement.

(a) Vary $a$, fix $\mathcal{R}_{\rm f}$ ($\mathcal{R}_{\rm f}=8$).  (b) Vary
$\mathcal{R}_{\rm f}$, fix $a$ ($a=0$).

Figure 2: Computation of
$S_{\rm{fd}}^{\infty}=\lim_{t\to\infty}S_{\rm{fd}}(t)$ via simulation of Eq.
10 with parameters $\mathcal{R}_{0}=2$, $p=0.25$, and $\delta=0.01$. The
initial conditions are $S=0.99$, $S_{\rm{fd}}=0$, and $I=0.01$. A bifurcation
occurs when $\kappa=1$ as defined by Eq. 30. When $\kappa<1$ the predicted
final size of the fearful susceptible group is $S_{\rm{fd}}^{\infty}=0$ while
for $\kappa>1$ it is given by Eq. 31. The results from the simulation are
plotted as a blue solid line, the predictions from the asymptotic methods are
plotted as a red dashed line and the predicted bifurcation point is a vertical
black dashed line.

We emphasize that the role of the long-term limit is to determine the
proportions of the susceptible population with and without fear that will
remain after the disease is no longer spreading through the population. The
reason fear can spread independently of the disease is because susceptible
people with fear are able to influence others to become afraid as well. The
role of fear is to adopt prophylactic behaviour (reduce contact through $p$),
however this parameter has an asymptotically small impact on disease burden
(see Fig. 3) after the disease outbreak. Therefore, for diseases with
parameters that follow the EDL, the use of fear to alter behaviour has very
little impact on disease dynamics and outcomes. However, the long-term limit
demonstrates that fear can remain in the system and thus people will adopt
prophylactic behaviours, but their efforts will be futile.

(a) Infected (b) Fearful susceptibles

Figure 3: Computation of infected individuals $I(t)$ via simulation of Eq. 10
with parameters $\mathcal{R}_{0}=2$, $\mathcal{R}_{\rm f}=8$, $a=0$, and
$\delta=0.01$ for different values of $p$. From these parameters $\kappa>1$ as
defined by Eq. 30 and thus the final susceptible group has a fearful
proportion. However, $p$ has little impact on the size of the infection. The
time scales on the two plots are different indicating that the dynamics of
fear settle on a much later timescale than those of the disease.

##  4 The Novel Disease Limit (NDL)

We demonstrated in Section 3 that the EDL cannot permit any behaviour-related
bifurcations in the infection outcomes of disease. We now consider
$\mathcal{R}_{\rm f}=\delta^{-1}\hat{\mathcal{R}}_{\rm f}$ and
$a=\delta^{-1}\hat{a}$ so that we are in the NDL. While numerically observing
discontinuity in epidemic prevalence, Epstein et al took parameters
$\beta_{\rm fd}\sim\beta$ which in the context of our scaling means
$\mathcal{R}_{\rm f}\gg 1$ and thus they were exploring what we now identify
as the NDL [14]. For simplicity and without loss of generality we will
consider $\hat{\mathcal{R}}_{\rm f}>\mathcal{R}_{0}-1$ so that following Eq.
18 the fearful susceptible population grows much faster than the disease. This
means that there will be a transient time when the fearless susceptible group
rapidly acquires fear of a novel disease before the epidemic grows
sufficiently. Therefore, for further simplicity we can consider the initial
conditions $S(0)=0$, $S_{\rm{fd}}(0)=1-\delta I_{0}$, and $I(0)=I_{0}$ noting
that, similar to the EDL in Section 3, a non-zero fearless susceptible
population does not impact the qualitative structure of the results. This
suggests that we write, $S=\delta x$, $S_{\rm{fd}}=1+\delta y$, and $I=\delta
z$ in Eq. 10 leading to the model for the NDL

(36a) $\displaystyle\dot{x}=$
$\displaystyle-\delta\mathcal{R}_{0}xz-\hat{\mathcal{R}}_{\rm
f}(1+\delta(y+z))x+1+\delta y-\hat{a}(x+y+z)(1+\delta y),$ (36b)
$\displaystyle\dot{y}=$ $\displaystyle-p\mathcal{R}_{0}(1+\delta
y)z+\hat{\mathcal{R}}_{\rm f}(1+\delta(y+z))x-(1+\delta
y)+\hat{a}(x+y+z)(1+\delta y),$ (36c) $\displaystyle\dot{z}=$
$\displaystyle\left(\mathcal{R}_{0}p-1+\delta\mathcal{R}_{0}(x+py)\right)z.$

Once again we pose an expansion

$x\sim x_{0}+\dots,\quad y\sim y_{0}+\dots,\quad z\sim z_{0}+\delta
z_{1}+\dots,$

where we note that the leading order behaviour in everything except the
disease variable will be sufficient. The leading order problem satisfies,

(37a) $\displaystyle\dot{x}_{0}=$ $\displaystyle-\hat{\mathcal{R}}_{\rm
f}x_{0}+1-\hat{a}(x_{0}+y_{0}+z_{0}),$ (37b) $\displaystyle\dot{y}_{0}=$
$\displaystyle-\mathcal{R}_{0}pz_{0}+\hat{\mathcal{R}}_{\rm
f}x-1+\hat{a}(x_{0}+y_{0}+z_{0}),$ (37c) $\displaystyle\dot{z}_{0}=$
$\displaystyle(\mathcal{R}_{0}p-1)z_{0},$

with initial conditions $x_{0}(0)=0$, $y_{0}(0)=-I_{0}$, and $z_{0}(0)=I_{0}$.
The leading order disease dynamics are decoupled from the fear transmission
similar to how the fearful susceptible group decoupled from disease dynamics
in the EDL of Section 3. This is why we are able to expand $z$ to an
additional order since the leading order of $x$ and $y$ will determine the
first correction to $z$. Thus,

(38) $\displaystyle z_{0}=I_{0}\textrm{e}^{{(\mathcal{R}_{\rm
p}-1)t}};\qquad\mathcal{R}_{\rm p}=\mathcal{R}_{0}p.$

Adding Eqs. 37a and 37b and integrating determines that

(39) $\displaystyle
x_{0}+y_{0}=-z_{0}-\hat{z}_{0};\qquad\hat{z}_{0}=\int_{0}^{t}z_{0}(s)\mathop{\mathrm{\mathstrut
d}}\nolimits\\!s=\frac{I_{0}}{\mathcal{R}_{\rm
p}-1}\left(\textrm{e}^{{(\mathcal{R}_{\rm p}-1)t}}-1\right),$

and this can be used to solve Eq. 37b yielding,

(40) $\displaystyle y(t)=y_{0}^{\infty}-A_{0}\textrm{e}^{{(\mathcal{R}_{\rm
p}-1)t}}-(y_{0}^{\infty}-A_{0}+I_{0})\textrm{e}^{{-\hat{\mathcal{R}}_{\rm
f}t}},$

where

(41) $\displaystyle A_{0}=\frac{\mathcal{R}_{\rm p}I_{0}}{\mathcal{R}_{\rm
p}-1}+\frac{\hat{a}I_{0}}{(\mathcal{R}_{\rm p}-1)(\mathcal{R}_{\rm
p}+\hat{\mathcal{R}}_{\rm f}-1)},\quad
y_{0}^{\infty}=\frac{(\hat{\mathcal{R}}_{\rm
f}+\hat{a})\frac{I_{0}}{\mathcal{R}_{\rm p}-1}-1}{\hat{\mathcal{R}}_{\rm f}}.$

We have defined $y_{0}^{\infty}$ as such since if $\mathcal{R}_{\rm p}<1$ then
$y_{0}\to y_{0}^{\infty}$ as $t\to\infty$. We note that $y_{0}^{\infty}<0$ is
fine since we have let $S_{\rm{fd}}=1+\delta y$ so the leading order solution
is already a correction from unity.

Since the disease is decoupled from fear transmission then the leading order
solutions determine the first correction for disease which satisfies,

(42) $\displaystyle\dot{z}_{1}=(\mathcal{R}_{\rm
p}-1)z_{1}-\mathcal{R}_{0}(z_{0}+\hat{z}_{0}+(1-p)y_{0})z_{0};\qquad
z_{1}(0)=0.$

This has solution,

(43) $\displaystyle
z_{1}(t)=-\mathcal{R}_{0}\left[U_{0}\textrm{e}^{{(\mathcal{R}_{\rm
p}-1)t}}+U_{1}\textrm{e}^{{-\hat{\mathcal{R}}_{\rm
f}t}}+u^{\infty}t-(U_{0}+U_{1})\right]z_{0},$

where

(44a) $\displaystyle U_{0}=$ $\displaystyle\frac{\mathcal{R}_{\rm
p}pI_{0}}{(\mathcal{R}_{\rm
p}-1)^{2}}-\frac{\hat{a}(1-p)I_{0}}{(\mathcal{R}_{\rm
p}-1)^{2}(\mathcal{R}_{\rm p}+\hat{\mathcal{R}}_{\rm f}-1)},$ (44b)
$\displaystyle U_{1}=$
$\displaystyle\left[\frac{\hat{a}I_{0}}{\hat{\mathcal{R}}_{\rm
f}^{2}(\mathcal{R}_{\rm p}+\hat{\mathcal{R}}_{\rm
f}-1)}-\frac{1}{\hat{\mathcal{R}}_{\rm f}^{2}}\right](1-p),$ (44c)
$\displaystyle u^{\infty}=$
$\displaystyle\frac{\hat{a}(1-p)I_{0}}{\hat{\mathcal{R}}_{\rm
f}(\mathcal{R}_{\rm p}-1)}-\frac{pI_{0}}{\mathcal{R}_{\rm
p}-1}-\frac{1-p}{\hat{\mathcal{R}}_{\rm f}}.$

We remark that in Eq. 38, yet another effective reproduction number,
$\mathcal{R}_{\rm p}$, emerges. The parameter $p$ is the prophylactic response
taken by the fearful susceptible class to reduce disease transmissibility.
Thus, if the disease is otherwise causing an outbreak ($\mathcal{R}_{0}>1$),
but the fearful behaviour response causes $\mathcal{R}_{\rm p}<1$ then the
prophylactic behaviour will mitigate disease impact. The leading order
response $z_{0}$ in Eq. 38 is an example of discontinuity in epidemic
prevalence as noted by Perra et. al in [32] since there will either be
exponential growth in the disease leading to a large proportion of disease
burden or exponential decay leading to minimal disease burden. Since this is a
leading order effect, we do not expect the true discontinuity in final disease
burden, but rather a rapid increase in the final recovered population past the
predicted bifurcation point $p_{c}=\mathcal{R}_{0}^{-1}$. We showcase this
bifurcation in Figure 4 by plotting $R_{\rm nat}^{\infty}$ from simulating Eq.
10 with $\hat{a}=0$, $\mathcal{R}_{\rm f}=10$, $\delta=0.01$, and varying $p$.
We take as initial conditions $S(0)=0$, $S_{\rm{fd}}(0)=0.99$ and $I(0)=0.01$.
We choose $\hat{\mathcal{R}}_{\rm f}$ such that $\hat{\mathcal{R}}_{\rm
f}>\mathcal{R}_{0}-1$ for all $\mathcal{R}_{0}$. The recovered group has an
inflection point very close to the predicted bifurcation value plotted in
dashed lines.

Figure 4: Computation of the final disease burden, $R_{\rm
nat}^{\infty}=\lim_{t\to\infty}R_{\rm nat}(t)$ from simulating Eq. 10 with
parameters $\hat{a}=0$, $\hat{\mathcal{R}}_{\rm f}=10$, and $\delta=0.01$
along with initial conditions $S(0)=0$, $S_{\rm{fd}}(0)=0.99$, $I(0)=0.01$.
The dashed lines are the predicted bifurcation points $\mathcal{R}_{\rm p}=1$
from Eq. 38.

The equilibrium analysis insures that the only disease steady state for the
infection is $I(t)=0$ and therefore unbounded exponential growth cannot
happen. The breakdown of the leading order infection Eq. 38 occurs at at time
when the exponential becomes large enough such that the assumption
$I\sim\mathcal{O}(\delta)$ fails. Therefore, when $\mathcal{R}_{\rm p}>1$, we
would rescale the equations to a long-time limit similar to Section 3.1 for
the EDL and analyze the final-size in disease burden. However, we instead
focus on the second, and more unintuitive, discontinuity in epidemic
prevalence associated to the self-limiting behaviour of the complacency
parameter $\hat{a}$. This is the unresolved mechanism noted by Perra et. al in
[32].

###  4.1 Weakly non-linear analysis near $\mathcal{R}_{0}^{-1}$

The complacency parameter $\hat{a}$ counters the prophylaxis of fearful
susceptibles since recovery from the disease demonstrates a reduction in risk
accelerating transfer away from being fearful. Thus, even if $\mathcal{R}_{\rm
p}<1$ it is possible that there exists a value of $\hat{a}$ such that the
disease burden does not remain small. We identify the role of the complacency
parameter on disrupting the disease burden final size near its natural
bifurcation point $p_{c}=\mathcal{R}_{0}^{-1}$ by letting

(45) $\displaystyle p=\mathcal{R}_{0}^{-1}-\sqrt{\delta}\hat{p},$

noting that we have taken a minus sign because we are primarily concerned with
the case $p<\mathcal{R}_{0}^{-1}$ ($\hat{p}>0$) since we have already
reconciled that when $p>\mathcal{R}_{0}^{-1}$ the exponential solutions Eq. 38
at leading order blow-up (see Fig. 4).

The expansion Eq. 45 is chosen as a distinguished limit of the system as
detailed in Appendix B. We also detail in that appendix how taking this
expansion requires a long-time analysis on a time scale $t\sim\delta^{-1/2}$
similar to the $\mathcal{O}(\delta)$ timescale needed for the fearful
susceptible population in the EDL as discussed in Section 3.1.

To analyze the long-term behaviour of the NDL we rescale time
$t=\delta^{-1/2}\tau$. Following Appendix B we must also rescale each of the
model variables $S=\sqrt{\delta}X$, $S_{\rm{fd}}=1+\sqrt{\delta}Y$, and
$I=\delta Z$ which transforms the model Eq. 10 to

(46a) $\displaystyle\sqrt{\delta}X_{\tau}=$
$\displaystyle-\hat{\mathcal{R}}_{\rm
f}X-\hat{a}(X+Y)+\sqrt{\delta}(1-\hat{\mathcal{R}}_{\rm
f}XY-\hat{a}(X+Y)Y-\hat{a}Z)+\mathcal{O}(\delta),$ (46b)
$\displaystyle\sqrt{\delta}Y_{\tau}=$ $\displaystyle\hat{\mathcal{R}}_{\rm
f}X+\hat{a}(X+Y)-\sqrt{\delta}(1-\hat{\mathcal{R}}_{\rm
f}XY-\hat{a}(X+Y)Y-\hat{a}Z+Z)+\mathcal{O}(\delta),$ (46c) $\displaystyle
Z_{\tau}=$
$\displaystyle(\mathcal{R}_{0}X+Y-\mathcal{R}_{0}\hat{p})Z+\mathcal{O}(\sqrt{\delta}),$

with initial conditions $X(0)=Y(0)=0$ and $Z(0)=I_{0}$.

Since $\delta\ll 1$ Eqs. 46a and 46b are both in a quasi-steady limit
yielding,

(47) $\displaystyle Y=-\left(\frac{\hat{\mathcal{R}}_{\rm
f}+\hat{a}}{\hat{a}}\right)X.$

However, this introduces a degeneracy into the system since both equations
yield the same result. This is resolved by adding Eqs. 46a and 46b,

(48) $\displaystyle X_{\tau}+Y_{\tau}=-Z+\mathcal{O(\sqrt{\delta})},$

which after combining with Eq. 47 and ignoring higher order terms yields,

(49) $\displaystyle\frac{\hat{\mathcal{R}}_{\rm f}}{\hat{a}}X_{\tau}=Z.$

Combining Eqs. 47, 49, and 46c allows us to determine an implicit relationship
between the fearless susceptible and infectious groups,

(50) $\displaystyle Z=\frac{\hat{\mathcal{R}}_{\rm
f}}{\hat{a}}\left(-\frac{\Gamma}{2}X^{2}-\mathcal{R}_{0}\hat{p}X+\frac{I_{0}\hat{a}}{\hat{\mathcal{R}}_{\rm
f}}\right);\qquad\Gamma=\frac{\hat{\mathcal{R}}_{\rm
f}-(\mathcal{R}_{0}-1)\hat{a}}{\hat{a}},$

where the initial conditions $X(0)=0$ and $Z(0)=I_{0}$ have been used. Further
combining Eq. 50 with Eq. 49 yields a first order non-linear differential
equation for $X$,

(51) $\displaystyle
X_{\tau}=-\frac{\Gamma}{2}X^{2}-\mathcal{R}_{0}\hat{p}X+\frac{I_{0}\hat{a}}{\hat{\mathcal{R}}_{\rm
f}},\quad X(0)=0.$

with solution,

(52) $\displaystyle
X=-\frac{\mathcal{R}_{0}\hat{p}}{\Gamma}+\frac{\sqrt{\theta}}{\hat{\mathcal{R}}_{\rm
f}\Gamma}\textrm{tanh}\left(\frac{\sqrt{\theta}}{2\hat{\mathcal{R}}_{\rm
f}}\tau+\hat{C}^{+}\right);\qquad\hat{C}^{+}=\textrm{arctanh}\left(\frac{\mathcal{R}_{0}\hat{\mathcal{R}}_{\rm
f}\hat{p}}{\sqrt{\theta}}\right),$

where

(53) $\displaystyle\theta=\mathcal{R}_{0}^{2}\hat{\mathcal{R}}_{\rm
f}^{2}\hat{p}^{2}+2I_{0}\hat{\mathcal{R}}_{\rm f}\hat{a}\Gamma.$

From Eq. 52 we have a classic final-size argument, namely that the final
susceptible group remains bounded,

$X^{\infty}=\lim_{\tau\to\infty}X=-\frac{\mathcal{R}_{0}\hat{p}}{\Gamma}+\frac{\sqrt{\theta}}{\hat{\mathcal{R}}_{\rm
f}\Gamma}$

and that the infection extinguishes. However, this is only true so long as
$\theta\geq 0$. Firstly, from Eq. 53 we note that if $\Gamma>0$ defined by Eq.
50 then necessarily $\theta>0$ always. Therefore, the solution can only begin
to break down when $\Gamma<0$ which occurs if $\hat{a}>\hat{a}_{c}$ where

(54) $\displaystyle\hat{a}_{c}=\frac{\hat{\mathcal{R}}_{\rm
f}}{\mathcal{R}_{0}-1}.$

However, even if $\hat{a}>\hat{a}_{c}$ then if $\hat{p}$ is sufficiently
large, $\theta$ can remain positive. Thus for $\theta$ to become negative and
for the solution Eq. 52 to break down we require $\hat{a}>\hat{a}_{c}$ and
$\hat{p}<\hat{p}_{c}$ where

(55)
$\displaystyle\hat{p}_{c}=\frac{1}{\mathcal{R}_{0}}\sqrt{\frac{2I_{0}\hat{a}|\Gamma|}{\hat{\mathcal{R}}_{\rm
f}}}.$

If $\hat{a}>\hat{a}_{c}$ and $\hat{p}<\hat{p}_{c}$ then the argument in the
hyperbolic tangent of Eq. 52 becomes complex and instead we can write the
solution as,

(56) $\displaystyle
X(\tau)=\frac{\mathcal{R}_{0}\hat{p}}{|\Gamma|}+\frac{\sqrt{|\theta|}}{\hat{\mathcal{R}}_{\rm
f}|\Gamma|}\tan\left(\frac{\sqrt{|\theta|}}{2\hat{\mathcal{R}}_{\rm
f}}\tau-\hat{C}^{-}\right),\qquad\hat{C}^{-}=\arctan\left(\frac{\mathcal{R}_{0}\hat{\mathcal{R}}_{\rm
f}\hat{p}}{\sqrt{|\theta|}}\right).$

This blows up when $\tau=\tau_{c}$ given by,

(57) $\displaystyle\tau_{c}=\frac{2\hat{\mathcal{R}}_{\rm
f}}{\sqrt{|\theta|}}\left(\frac{\pi}{2}-\hat{C}^{-}\right),$

at which point a new rescaling would be needed to resolve the relaxation to
equilibrium.

We demonstrate the bifurcation behaviour around $\hat{p}_{c}$ in Fig. 5 where
we simulate Eq. 10 with initial conditions $S(0)=0$, $S_{\rm{fd}}(0)=0.99$ and
$I(0)=0.01$ while taking $\delta=0.01$. For Figs. 5a and 5b we fix
$\mathcal{R}_{0}=2$, $\hat{a}=2$, and $p=0.48$ ($\hat{p}=0.2$) and vary
$\hat{\mathcal{R}}_{\rm f}$ taking $\hat{\mathcal{R}}_{\rm f}=3,4,5,6$. This
results in $\Gamma=0.5,1,1.5$, and $2$ respectively which are all positive and
thus $\theta>0$ defined by Eq. 53. Therefore, the solutions for the fearless
susceptibles should follow $S=\sqrt{\delta}X$ given by Eq. 52 while the
fearful susceptibles are given by $S_{1}=1+\sqrt{\delta}Y$ given by Eq. 47. We
note excellent agreement between the simulation (solid lines) and analytical
approximations (dashed lines) improving with increasing
$\hat{\mathcal{R}}_{\rm f}$. This improvement is unsurprising since, for
example, when $\hat{\mathcal{R}}_{\rm f}=3$ then $\Gamma=0.5\sim\sqrt{\delta}$
and so other asymptotic structures are involved.

For Figs. 5c and 5d we fix $\mathcal{R}_{0}=2$, $\hat{\mathcal{R}}_{\rm f}=3$,
$\hat{a}=8$, and $p=0.44$ ($\hat{p}=0.6$). With these values
$\Gamma=-0.625<0$. Furthermore from Eq. 55 $\hat{p}_{c}=0.9129$ and since
$\hat{p}<\hat{p}_{c}$ then $\theta<0$ and so we expect a bifurcation.
Therefore, the solutions for the fearless susceptibles now should follow
$S=\sqrt{\delta}X$ given now by Eq. 56. Indeed we observe this in Fig. 5c
where the agreement is excellent until approximately $\tau_{c}$ (black dashed
line) where there is a change in the solution behaviour and the asymptotic
approximation fails. We show how this leads to an epidemic in the infected
class in Fig. 5d where $I=\delta Z$ is given by Eq. 50.

(a) $\theta>0$, $p=0.48$, vary $\hat{\mathcal{R}}_{\rm f}$. (b) $\theta>0$,
$p=0.48$, vary $\hat{\mathcal{R}}_{\rm f}$.

(c) $\theta<0$, $p=0.44$, $\hat{\mathcal{R}}_{\rm f}=3$.  (d) $\theta<0$,
$p=0.44$, $\hat{\mathcal{R}}_{\rm f}=3$.

Figure 5: Comparison of simulations of Eq. 10 with parameters $\hat{a}=2$,
$\mathcal{R}_{0}=2$, and $\delta=0.01$ along with initial conditions $S(0)=0$,
$S_{\rm{fd}}(0)=0.99$, $I(0)=0.01$. The dashed curves are the asymptotic
approximation valid on a time scale $t\sim\delta^{-1/2}$ given by Eq. 52 (Fig.
5a), Eq. 47 (Fig. 5b), Eq. 56 (Fig. 5c) and Eq. 50 (Fig. 5d). The vertical
black dashed line indicates the approximation $\tau_{c}$ given by Eq. 57 for
the time when the asymptotic solution fails if $\theta<0$.

Overall, we have deciphered the mechanism first reported by Perra et al in
[32] concerning discontinuity in epidemic prevalence. When the prophylactic
behaviour of the fearful susceptible population is sufficient to prevent an
epidemic ($\mathcal{R}_{\rm p}<1$) then the self-limiting behaviour of
complacency causes an outbreak to emerge when $\hat{a}>\hat{a}_{c}$ and
$p>\mathcal{R}_{0}^{-1}-\sqrt{\delta}\hat{p}_{c}$ with $\hat{p}_{c}$ given by
Eq. 55 as demonstrated in Fig. 5d. We demonstrate the discontinuity in
epidemic prevalence in Fig. 6 simulating Eq. 10 with $\mathcal{R}_{0}=2$,
$\hat{\mathcal{R}}_{\rm f}=3$, and $\delta=0.01$ while varying $\hat{a}$. We
plot the predicted bifurcation point
$\mathcal{R}_{0}^{-1}-\sqrt{\delta}\hat{p}_{c}$ from Eq. 55 as vertical dashed
lines. We see that the approximations worsen as $\hat{a}$ increases which is
expected since, for example if $\hat{a}=10=1/\sqrt{\delta}$ then other
asymptotic considerations need to be made.

Figure 6: The final epidemic prevalence $R^{\infty}$ simulated from Eq. 10
with parameters $\mathcal{R}_{0}=2$, $\hat{\mathcal{R}}_{\rm f}=3$ and
$\delta=0.01$ along with initial conditions $S(0)=0$, $S_{\rm{fd}}(0)=0.99$,
$I(0)=0.01$ while varying $\hat{a}$ and $p$. The dashed curves are the
predicted points of discontinuity in epidemic prevalence from Eq. 55.

We plot the bifurcation curve $p_{c}(\hat{a})$ numerically generated from
simulating Eq. 10 for $\mathcal{R}_{0}=2$ and $\hat{\mathcal{R}}_{\rm f}=3$ to
that predicted from Eq. 55 in Fig. 7a (for $\delta=0.01$) and Fig. 7b (for
$\delta=0.001$) showing excellent agreement within $\mathcal{O}(\delta)$, the
order considered here.

(a) $\delta=0.01$.  (b) $\delta=0.001$.

Figure 7: The epidemic prevalence discontinuity bifurcation curve simulated
(solid lines) from Eq. 10 with parameters $\mathcal{R}_{0}=2$ and
$\hat{\mathcal{R}}_{\rm f}=3$ along with initial conditions $S(0)=0$,
$S_{\rm{fd}}(0)=0.99$, $I(0)=0.01$. The dashed lines are the predicted
bifurcation curves given by
$p_{c}(\hat{a})=\mathcal{R}_{0}^{-1}-\sqrt{\delta}\hat{p}_{c}$ with
$\hat{p}_{c}$ given by Eq. 55.

Finally, we plot the infected class $I(t)$ varying $p$ in Fig. 8. We see that
when $p<p_{c}$ then there is no outbreak as predicted (blue curves). However,
when $p>p_{c}$ not only does the disease outbreak occur, it is accompanied by
a second outbreak (red curves). We also plot in black the infection curve if
behaviour is not included in the model at all (no $S_{\rm{fd}}$ compartment).
From Fig. 8 we see that the role of prophylactic behaviour when an outbreak
occurs is that the initial wave has a reduced burden compared to when there is
no prophylaxis. However, as this behaviour wanes then additional waves occur.
We leave the analysis of additional waves to future work.

Figure 8: Simulations of $I(t)$ from Eq. 10 with parameters
$\mathcal{R}_{0}=2$, $\hat{\mathcal{R}}_{\rm f}=3$, $\hat{a}=8$, and
$\delta=0.01$ along with initial conditions $S(0)=0$, $S_{\rm{fd}}(0)=0.99$,
$I(0)=0.01$. The blue dashed curves are trajectories for $p<p_{c}$ the
bifurcation curve in Fig. 5 while the red solid curves are for $p>p_{c}$. The
black curve is the simulation of a basic SIR model without behavioural
intervention (i.e. the fearful susceptible compartment $S_{\rm{fd}}$ is not
present).

##  5 Discussion and Conclusions

We have considered a behaviour-modified susceptible-infected-removed
mathematical model within which disease-susceptible individuals can acquire
fear of the disease and therefore reduce their contact behaviour and hence
their susceptibility. Such behaviour considerations are important since the
traditional SIR model with constant parameters can only permit a single
outbreak, however recurrent waves of disease during a single outbreak are
often realized [29]. We observed that the combination of prophylactic and
complacent behaviours can lead to multiple disease outbreaks in Fig. 8.

Often the effect of behaviour on disease transmission is captured through
estimation of a time-dependent effective reproduction number [20]. However,
estimates of this can vary dependent on the assumptions, data, and statistical
methods used to calculated it [8]. Thus, more mechanistic considerations of
behaviour can identify predictive ways in which behaviour affects disease
transmission through meaningful, identifiable, and measurable parameters as
well as a way to inform parameters by incorporating behaviour and other
sociological data.

Beginning with a naive and fearless susceptible population introduced to a
small initial infection, we showed that two distinct dynamic limits emerge.
The first one occurs when the reproduction numbers of fear and disease are
comparable meaning that fear is acquired on a timescale much slower than the
infection time. We called such a limit the Establish Disease Limit because
there is an implication of familiarity with the disease and its consequences
since the behaviour response trails increasing case numbers. The second limit
was one in which the behaviour response was rapid and preceded growth in case
numbers. Correspondingly the fear reproduction number is much larger than that
of the disease implying that there is a relatively unknown understanding of
the disease burden. As such we called this the Novel Disease Limit.

The explicit inclusion of behaviour in the model allowing for two distinct
dynamic limits is very insightful. It shows how important human behaviour is
in assessing disease impacts and burden reduction strategy. It also provides
an early-warning system for policy makers since data from the onset of an
outbreak can be used to estimate model parameters. These parameter values
would indicate if the population response was in the Established or Novel
disease limits providing insight into mitigation strategies that would be
effective.

The EDL analyzed in Section 3 showed that a disease outbreak will occur if the
disease basic reproduction number $\mathcal{R}_{0}$ is greater than one. This
is an unsurprising result because this limit has a delayed behavioural
response. Once the epidemic ends then there is a long-time limit on the
timescale of fear recovery $\gamma_{\rm f}^{-1}$ where fear response becomes
important. At this point, since the disease spread has ended, the model
behaves as a susceptible-infected-susceptible (SIS) model with fear as the
“infection”.

A characteristic result of the SIS model is that the disease either vanishes
or becomes endemic [6]. We confirmed this feature in Fig. 2 where we showed
that the presence of fear at steady state was dependent on the fear
reproduction number $\mathcal{R}_{\rm f}$ and the fear-loss due to
complacency, $a$. The endemicity of infection (analogous to fear in our case)
implies that a weak fear response can result in people retaining fear even
after the outbreak is over. Such an endemic fear state has been observed for
example during COVID-19 when people continued wearing masks while the rate of
disease spread was slow [25].

Since the behavioural response to the disease is delayed in the EDL then the
population is not taking prophylactic action when needed at the onset of the
outbreak, but are taking it when the disease outbreak is over and it is not
needed. We showed in Fig. 3a that because of this inverted response,
prophylactic behaviour has limited ability to mitigate epidemic prevalence in
the EDL. Having a fearful susceptible population will reduce any resurgent
outbreaks since they will be adopting prophylactic behaviour. However, this
may only have value for a short duration between sequential outbreaks as
natural and complacency fear removal mechanisms will reduce the fearful
population. There may also be negative consequences with retaining fear in the
population. For example, people may observe that their behavioural changes are
unimpactful and resist behavioural change in future outbreaks. There may also
be civil unrest, particularly if the behavioural change is imposed through
government action [10].

Through a weakly nonlinear analysis in the NDL presented in Section 4 we were
able to resolve the bifurcation observed by Perra et al. in [32] that even if
the prophylactic strength is such that $\mathcal{R}_{\rm p}<1$ then it may not
be enough to suppress an outbreak. This is possible because, while the role of
$p$ is to reduce disease transmission, there is a second parameter, $\hat{a}$
which is a behaviour complacency rate. This causes people to relax their
behaviour to normal conditions and hence in the context of the model, lose
fear. This rate of fear loss is proportional to the number of people who have
recovered, thus as the burden and consequences of the disease become apparent
and people begin to get better, there is less urgency to be afraid and people
relax. Overall then, there is a competition between mitigating behaviour that
reduces outbreaks and complacent behaviour which encourages them. This leads
to the discontinuity in epidemic prevalence observed in Fig. 6.

Epidemiologically, the discontinuity in epidemic prevalence occurs because
when complacency is sufficiently high then a significant proportion of people
stop contributing to the prophylactic behaviour. Those that remain must
increase the strength of their prophylaxis to overcome this decrease in
population. Thus, since $1-p$ represents a contact reduction, $p$ must
decrease to continue to prevent an epidemic.

Understanding the tipping point in epidemic prevalence is of great value to
public health. If there is an estimate of population prophylaxis and sentiment
to adherence behaviour or complacency then our analysis can predict how well
the disease is being managed and how close we are to an epidemic outbreak.
Understanding the mechanism can also help with the development of public
health messaging. If the mitigating behaviour is sufficient such that the
reproduction number falls below one then when disease burden continues to
climb, people may feel that their efforts have no value and fatigue sets in
[34]. However, our result allows a quantification of the required prophylactic
behaviour to suppress disease spread which could encourage adherence.

Having identified two distinct disease limits, it is of interesting future
work to consider the transition between novel and established diseases such as
when a new disease becomes endemic in the population. In fact this transition
in human behaviour could lead to a robust definition of an endemic disease.
Finally, we remark that in the model analyzed here, acquisition of fear is
directly proportional to the infected class which mimics the mass-action
kinetics between susceptible and infected individuals in a standard SIR model.
While that assumption can be justified epidemiologically because two
individuals in close approximation can spread a disease unaware of the health
status of the other, fear requires a cognizant awareness. Thus, the model
assumes that if infection is driving fear then either susceptible people know
the infection status of those they come in contact with or are being provided
accurate real-time information of the number of infections. In reality it is
likely that neither of things are true. Therefore, it will be interesting to
explore other pathways of fear acquisition including under-reporting of cases
or intentional misinformation campaigns.

##  Appendix A Solution Breakdown in the Established Disease Limit

For the EDL we solve Eq. 19 by expanding in an asymptotic expansion for
$\delta$,

$S\sim x_{0}+\delta x_{1},\quad S_{\rm{fd}}\sim\delta
y_{0}+\delta^{2}y_{1},\quad I\sim z_{0}+\delta z_{1}.$

The leading order solutions are given by Eq. 22 and Eq. 23 and as $t\to\infty$
have values $z_{0}^{\infty}=0$, $x_{0}^{\infty}$ given by Eq. 24 and
$y_{0}^{\infty}$ given by Eq. 23. The next order equations are Eq. 25 which if
we consider a long-time where $x_{0}\to x_{0}^{\infty}$, $y_{0}\to
y_{0}^{\infty}$, and $z_{0}\to 0$ then these differential equations are
approximately given by

(58) $\displaystyle\dot{x}_{1}\approx$
$\displaystyle-\mathcal{R}_{0}x_{0}^{\infty}z_{1},$ (59)
$\displaystyle\dot{y}_{1}\approx$ $\displaystyle(\mathcal{R}_{\rm
f}x_{0}^{\infty}-\mathcal{R}_{0}py_{0}^{\infty})z_{1}+\mathcal{R}_{\rm
f}x_{0}^{\infty}y_{0}^{\infty}-y_{0}^{\infty}-a(1-x_{0}^{\infty})y_{0}^{\infty},$
(60) $\displaystyle\dot{z}_{1}\approx$
$\displaystyle(\mathcal{R}_{0}x_{0}^{\infty}-1)z_{1}.$

The last of these has solution,

(61) $\displaystyle
z_{1}=C_{z}\textrm{e}^{{-(1-\mathcal{R}_{0}x_{0}^{\infty})t}},$

for constant $C_{z}$ which would be determined from early-time considerations.
Since $x_{0}^{\infty}<\mathcal{R}_{0}^{-1}$ then this decays and we can assume
in the long-time limit that the disease has gone extinct. We can then solve
Eq. 59 yielding,

(62) $\displaystyle
y_{1}=C_{y}+\frac{(\mathcal{R}_{0}py_{0}^{\infty}-\mathcal{R}_{\rm
f}x_{0}^{\infty})}{1-\mathcal{R}_{0}x_{0}^{\infty}}C_{z}\textrm{e}^{{-(1-\mathcal{R}_{0}x_{0}^{\infty})t}}+y_{0}^{\infty}\left(\mathcal{R}_{\rm
f}x_{0}^{\infty}-1-1(1-x_{0}^{\infty})\right)t,$

for another constant $C_{y}$. Problematically, we note that Eq. 62 grows
without bound due to the linear term. Since the expansion is $y_{0}+\delta
y_{1}$ then this sequence of terms will no longer be asymptotic when
$t\sim\mathcal{O}(\delta^{-1})$. Therefore we need to consider Eq. 10 with a
new timescale $t=\delta^{-1}\tau$ where $I=0$ to determine the solution
behaviour of $S$ and $S_{\rm{fd}}$ in a long-time limit. As $\tau\to 0$ these
new solutions must match to the solutions for $t\sim\mathcal{O}(1)$ as
$t\to\infty$ (cf. [21]). Therefore,

(63) $\displaystyle\lim_{t\to\infty}S(t)=$ $\displaystyle
x_{0}^{\infty}+\delta x_{1}^{\infty}\sim\lim_{\tau\to 0}S(\tau)$ (64)
$\displaystyle\lim_{t\to\infty}S_{\rm{fd}}(t)=$ $\displaystyle\delta
y_{0}^{\infty}\sim\lim_{\tau\to 0}S_{\rm{fd}}(\tau).$

Since these are just constant then this matching is equivalent to taking
$S=x_{0}^{\infty}+\delta x_{1}^{\infty}$ and $S_{\rm{fd}}=\delta
y_{0}^{\infty}$ when $\tau=0$ which are the initial conditions used in Section
3.1.

##  Appendix B Solution Breakdown in the Novel Disease Limit

The leading order behaviour in the Novel Disease Limit as discussed in Section
4 has solutions with an exponential argument $(\mathcal{R}_{0}p-1)$ and thus
the solution behaves very differently around a critical point
$p=\mathcal{R}_{0}^{-1}$. To see how this appears more directly in the model,
consider starting with a small initial infection $I\sim\mathcal{O}(\delta)$
and looking at a long-time limit to Eq. 10. The first two equations relate the
susceptible compartments which produces two possible quasi-steady states,

(65) $\displaystyle S_{\rm{fd}}=0;\qquad
S_{\rm{fd}}=1-\frac{\hat{\mathcal{R}}_{\rm f}+\hat{a}}{\hat{\mathcal{R}}_{\rm
f}}S$

and we will focus on the latter fear-present case since the leading order
fearful behaviour in Section 4 given by Eq. 40 is non-zero. The quasi-steady
equation for infection Eq. 10c is,

(66) $\displaystyle 0=f(S,S_{\rm{fd}})I;\qquad
f(S,S_{\rm{fd}})=\mathcal{R}_{0}(S+pS_{\rm{fd}})-1.$

From Eq. 66 the only quasi-steady solution is $I=0$ unless
$f(S,S_{\rm{fd}})=0$ which occurs when

(67) $\displaystyle S_{\rm{fd}}=\frac{1-\mathcal{R}_{0}S}{\mathcal{R}_{0}p}.$

This indicates a transition in the dynamics of infection from decay to growth.
Combining Eqs. 65 and 67 then the dynamics of the infection will begin to
transition when $S=S^{\dagger}$ given by

(68) $\displaystyle S^{\dagger}=\frac{(\mathcal{R}_{\rm
p}-1)\hat{a}}{\mathcal{R}_{0}((\hat{\mathcal{R}}_{\rm f}+\hat{a})p-\hat{a}}.$

Initially in Section 4 we consider, without loss of generality, that
$S_{\rm{fd}}\approx 1$ and $S\approx 0$ and so
$f(S,S_{\rm{fd}})\approx\mathcal{R}_{0}p-1$. Assuming that
$\mathcal{R}_{0}p<1$ then the long term behaviour of fearful susceptibles is
$S_{\rm{fd}}\sim 1+\delta y_{0}^{\infty}$ given by Eq. 41 while the long term
behaviour of fearless susceptibles is $S\sim\delta x_{0}^{\infty}$ where

(69) $\displaystyle x_{0}^{\infty}=\frac{I_{0}\hat{a}+1-\mathcal{R}_{\rm
p}}{\hat{\mathcal{R}}_{\rm f}(1-\mathcal{R}_{\rm p})}.$

Since the infection is initially decaying then we are interested in the values
of $p$ where $S=\delta x_{0}^{\infty}=S^{\dagger}$ given by Eq. 68 since that
is when the infection behaviour can transition to growth. Equating these leads
to a rational expression with quadratic numerator,

(70)
$\displaystyle\frac{a_{2}p^{2}-a_{1}p+a_{0}}{\mathcal{R}_{0}\hat{\mathcal{R}}_{\rm
f}(\mathcal{R}_{\rm p}-1)(\mathcal{R}_{\rm f}p-\hat{a}(1-p))}=0,$

where

(71a) $\displaystyle a_{0}=$ $\displaystyle\hat{a}(\hat{\mathcal{R}}_{\rm
f}-\delta\mathcal{R}_{0}(1+I_{0}\hat{a})),$ (71b) $\displaystyle a_{1}=$
$\displaystyle\mathcal{R}_{0}(2\hat{a}\hat{\mathcal{R}}_{\rm
f}-\delta(\hat{\mathcal{R}}_{\rm f}+\hat{a}(I_{0}\hat{\mathcal{R}}_{\rm
f}+\mathcal{R}_{0}+1)+I_{0}\hat{a}^{2}),$ (71c) $\displaystyle a_{2}=$
$\displaystyle\mathcal{R}_{0}^{2}(\hat{\mathcal{R}}_{\rm
f}\hat{a}-\delta(\hat{\mathcal{R}}_{\rm f}+\hat{a}).$

When $\delta>0$ then Eq. 70 has two distinct solutions for $p$. However, when
$\delta=0$ then the quadratic has a single degenerate root
$p=\mathcal{R}_{0}^{-1}$. This confirms that $\mathcal{R}_{\rm p}\approx 1$
triggers the bifurcation, however the degenerate nature suggests the
corrections are expanded in powers of $\sqrt{\delta}$ [21]. Thus when
analyzing the behaviour near $p=\mathcal{R}_{0}^{-1}$ we consider this
distinguished limit and let,

(72) $\displaystyle p=\mathcal{R}_{0}^{-1}-\sqrt{\delta}\hat{p}.$

We substitute Eq. 72 in the leading order behaviour $S_{\rm{fd}}=1+\delta
y_{0}$ and $I=\delta z_{0}$ with $y_{0}$ and $z_{0}$ given by Eqs. 40 and 38
respectively and expand for $\delta\ll 1$ yielding,

(73a) $\displaystyle y_{0}\sim$ $\displaystyle-\frac{\hat{\mathcal{R}}_{\rm
f}+\hat{a}}{\hat{\mathcal{R}}_{\rm
f}}I_{0}t+\frac{I_{0}\hat{a}}{\hat{\mathcal{R}}_{\rm
f}^{2}}-\frac{1}{\hat{\mathcal{R}}_{\rm
f}}-I_{0}-\left(\frac{I_{0}\hat{a}}{\hat{\mathcal{R}}_{\rm
f}^{2}}-\frac{1}{\hat{\mathcal{R}}_{\rm
f}}\right)\textrm{e}^{{-\hat{\mathcal{R}}_{\rm f}t}},$ (73b) $\displaystyle
z_{0}\sim$ $\displaystyle I_{0}(1-\mathcal{R}_{0}\sqrt{\delta}\hat{p}t).$

Similarly from the relation Eq. 39 then

(74) $\displaystyle
x_{0}=-y_{0}-z_{0}-\hat{z}_{0},\qquad\hat{z}_{0}\underset{\delta\ll
1}{\approx}I_{0}\left(t-\sqrt{\delta}\frac{\mathcal{R}_{0}\hat{p}}{2}t^{2}\right).$

We see in Eq. 73b that the asymptotic sequence breaks down if
$t\sim\mathcal{O}(\delta^{-1/2})$ and so we let $t=\delta^{-1/2}\tau$ which we
substitute into Eqs. 73 and 74 producing to leading order,

(75a) $\displaystyle S\sim\delta x_{0}\approx$
$\displaystyle\delta^{1/2}\hat{a}I_{0}\tau,$ (75b) $\displaystyle
S_{\rm{fd}}\sim 1+\delta y_{0}\approx$ $\displaystyle
1-\delta^{1/2}\frac{\hat{\mathcal{R}}_{\rm f}+\hat{a}}{\hat{\mathcal{R}}_{\rm
f}}I_{0}\tau,$ (75c) $\displaystyle I\sim\delta z_{0}\approx$
$\displaystyle\delta I_{0}(1-\mathcal{R}_{0}\hat{p}\tau),$

which suggests that for this long-time limit we expand,

(76) $\displaystyle S\sim\delta^{1/2}X,\quad S_{\rm{fd}}\sim
1+\delta^{1/2}Y,\quad I\sim\delta Z.$

This is precisely the scaling used in Section 4.1 to perform the weakly non-
linear analysis in the long-time limit near $p=\mathcal{R}_{0}^{-1}$. Similar
to the EDL in Appendix A, matching will necessitate that as $\tau\to 0$ then
$X\sim Y\sim 0$ and $Z\sim I_{0}$.

## References

  * [1] C. Bauch, A. D’Onofrio, and P. Manfredi, Behavioral epidemiology of infectious diseases: An overview, Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases, (2013), pp. 1–19, https://doi.org/10.1007/978-1-4614-5474-8_1.
  * [2] B. Baumgaertner, J. E. Carlisle, and F. Justwan, The influence of political ideology and trust on willingness to vaccinate, PLoS ONE, 13 (2018), https://doi.org/10.1371/JOURNAL.PONE.0191728.
  * [3] B. Baumgaertner, B. J. Ridenhour, F. Justwan, J. E. Carlisle, and C. R. Miller, Risk of disease and willingness to vaccinate in the united states: A population-based survey, PLoS Medicine, 17 (2020), https://doi.org/10.1371/JOURNAL.PMED.1003354.
  * [4] B. O. Baumgaertner, P. A. Fetros, S. M. Krone, and R. C. Tyson, Spatial opinion dynamics and the effects of two types of mixing, Physical Review E, 98 (2018), https://doi.org/10.1103/PhysRevE.98.022310.
  * [5] B. O. Baumgaertner, R. C. Tyson, and S. M. Krone, Opinion strength influences the spatial dynamics of opinion formation, Journal of Mathematical Sociology, 40 (2016), pp. 207–218, https://doi.org/10.1080/0022250X.2016.1205049.
  * [6] F. Brauer, C. Castillo-Chavez, and Z. Feng, Mathematical Models in Epidemiology, vol. 69, Springer New York, 2019, https://doi.org/10.1007/978-1-4939-9828-9, https://link.springer.com/10.1007/978-1-4939-9828-9.
  * [7] F. Brauer, P. Van den Driessche, J. Wu, and L. J. Allen, Mathematical epidemiology, vol. 1945, Springer, 2008.
  * [8] E. K. Brockhaus, D. Wolffram, T. Stadler, M. Osthege, T. Mitra, J. M. Littek, E. Krymova, A. J. Klesen, J. S. Huisman, S. Heyder, et al., Why are different estimates of the effective reproductive number so different? a case study on covid-19 in germany, PLOS Computational Biology, 19 (2023), p. e1011653.
  * [9] J. K. Carpenter, M. Pinaire, and S. G. Hofmann, From extinction learning to anxiety treatment: Mind the gap, Brain sciences, 9 (2019), p. 164.
  * [10] R. Cordell, R. M. Wood, and T. M. Wright, Disease and dissent: Epidemics as a catalyst for social unrest, Global Studies Quarterly, 3 (2023), p. ksad031.
  * [11] P. V. den Driessche, Reproduction numbers of infectious disease models, Elsevier, 2 (2017), pp. 282–303, https://www.sciencedirect.com/science/article/pii/S2468042717300209.
  * [12] P. V. D. Driessche and J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical Biosciences, 180 (2002), pp. 29–48, https://doi.org/10.1016/S0025-5564(02)00108-6.
  * [13] D. P. Durham and E. A. Casman, Incorporating individual health-protective decisions into disease transmission models: a mathematical framework, Journal of the Royal Society Interface, 9 (2012), pp. 562–570, https://royalsocietypublishing.org/doi/abs/10.1098/rsif.2011.0325.
  * [14] J. M. Epstein, E. Hatna, and J. Crodelle, Triple contagion: A two-fears epidemic model, Journal of the Royal Society Interface, 18 (2021), https://doi.org/10.1098/RSIF.2021.0186.
  * [15] J. M. Epstein, J. Parker, D. Cummings, and R. A. Hammond, Coupled contagion dynamics of fear and disease: Mathematical and computational explorations, PLoS ONE, 3 (2008), https://doi.org/10.1371/JOURNAL.PONE.0003955.
  * [16] S. Funk, M. Salathé, and V. A. Jansen, Modelling the influence of human behaviour on the spread of infectious diseases: A review, Journal of the Royal Society Interface, 7 (2010), pp. 1247–1256, https://doi.org/10.1098/RSIF.2010.0142, https://drive.google.com/file/u/1/d/1CfQ4HYVPN656gTI_U7q6S_JR_ggM6RgP/view?usp=embed_facebook.
  * [17] G. P. Garnett, The geographical and temporal evolution of sexually transmitted disease epidemics, Sexually transmitted infections, 78 (2002), pp. i14–i19.
  * [18] G. P. Garnett, Role of herd immunity in determining the effect of vaccines against sexually transmitted disease, The Journal of infectious diseases, 191 (2005), pp. S97–S106.
  * [19] A. Glaubitz, F. F. P. of the Royal Society A, and undefined 2020, Oscillatory dynamics in the dilemma of social distancing, Proceedings of the Royal Society A, 476 (2020), https://doi.org/10.1098/rspa.2020.0686, https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2020.0686.
  * [20] K. M. Gostic, L. McGough, E. B. Baskerville, S. Abbott, K. Joshi, C. Tedijanto, R. Kahn, R. Niehus, J. A. Hay, P. M. D. Salazar, J. Hellewell, S. Meakin, J. D. Munday, N. I. Bosse, K. Sherrat, R. N. Thompson, L. F. White, J. S. Huisman, J. Scire, S. Bonhoeffer, T. Stadler, J. Wallinga, S. Funk, M. Lipsitch, and S. Cobey, Practical considerations for measuring the effective reproductive number, rt, PLoS Computational Biology, 16 (2020), https://doi.org/10.1371/JOURNAL.PCBI.1008409.
  * [21] E. J. Hinch, Perturbation Methods, Cambridge University Press, 10 1991, https://doi.org/10.1017/CBO9781139172189, https://www.cambridge.org/core/books/perturbation-methods/78E3D7607E441E4BCAE84698DE91D3BC.
  * [22] Y. Iwasa and R. Hayashi, Waves of infection emerging from coupled social and epidemiological dynamics, Journal of Theoretical Biology, 558 (2023), https://www.sciencedirect.com/science/article/pii/S0022519322003575.
  * [23] W. O. Kermack and A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 115 (1927), pp. 700–721, https://doi.org/10.1098/RSPA.1927.0118.
  * [24] Y. Kim, C. W. Won, S. Kim, B. Kim, M. Kim, E. Jeong, J. Yang, and H. Lee, Perceived recovery time from common cold as a possible indicator of physical resilience, Annals of geriatric medicine and research, 25 (2021), p. 204.
  * [25] M. Kwon and W. Yang, Mask-wearing behaviors after two years of wearing masks due to covid-19 in korea: A cross-sectional study, International Journal of Environmental Research and Public Health, 19 (2022), https://doi.org/10.3390/IJERPH192214940, /pmc/articles/PMC9691200//pmc/articles/PMC9691200/?report=abstracthttps://www.ncbi.nlm.nih.gov/pmc/articles/PMC9691200/.
  * [26] B. Liu, D. Jayasundara, V. Pye, T. Dobbins, G. J. Dore, G. Matthews, J. Kaldor, and P. Spokes, Whole of population-based cohort study of recovery time from covid-19 in new south wales australia, The Lancet Regional Health–Western Pacific, 12 (2021).
  * [27] J. Molla, S. Farhang-Sardroodi, I. R. Moyles, and J. M. Heffernan, Pharmaceutical and non-pharmaceutical interventions for controlling the covid-19 pandemic, Royal Society Open Science, (2023), https://doi.org/10.1098/rsos.230621.
  * [28] I. R. Moyles, J. M. Heffernan, and J. D. Kong, Cost and social distancing dynamics in a mathematical model of covid-19 with application to ontario, canada, Royal Society open science, 8 (2021), https://doi.org/10.1098/rsos.201770, https://royalsocietypublishing.org/doi/abs/10.1098/rsos.201770.
  * [29] A. Mummert and H. Weiss, Get the news out loudly and quickly: The influence of the media on limiting emerging infectious disease outbreaks, PLOS ONE, 8 (2013), p. e71692, https://doi.org/10.1371/JOURNAL.PONE.0071692, https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0071692.
  * [30] K. Myers, K. Ressler, and M. Davis, Different mechanisms of fear extinction dependent on length of time since fear acquisition, Learning & Memory, 13 (2006), pp. 216–223, https://doi.org/10.1101/lm.119806, http://learnmem.cshlp.org/content/13/2/216.short.
  * [31] L. G. Nardin, C. R. Miller, B. J. Ridenhour, S. M. Krone, P. Joyce, and B. O. Baumgaertner, Planning horizon affects prophylactic decision-making and epidemic dynamics, PeerJ, 4 (2016), https://peerj.com/articles/2678/.
  * [32] N. Perra, D. Balcan, B. Gonçalves, and A. Vespignani, Towards a characterization of behavior-disease models, PLoS ONE, 6 (2011), https://doi.org/10.1371/JOURNAL.PONE.0023084.
  * [33] N. Perra and A. Vespignani, Modeling contact and mobility based social response to the spreading of infectious diseases, Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases, (2013), pp. 103–123, https://doi.org/10.1007/978-1-4614-5474-8_7.
  * [34] A. Petherick, R. Goldszmidt, E. B. Andrade, R. Furst, T. Hale, A. Pott, and A. Wood, A worldwide assessment of changes in adherence to covid-19 protective behaviours and hypothesized pandemic fatigue, Nature Human Behaviour, 5 (2021), pp. 1145–1160.
  * [35] P. Poletti, M. Ajelli, and S. Merler, The effect of risk perception on the 2009 h1n1 pandemic influenza dynamics, PLoS ONE, 6 (2011), https://doi.org/10.1371/JOURNAL.PONE.0016460.
  * [36] F. Raeder, C. J. Merz, J. Margraf, and A. Zlomuzica, The association between fear extinction, the ability to accomplish exposure and exposure therapy outcome in specific phobia, Scientific reports, 10 (2020), p. 4288.
  * [37] R. C. Tyson, N. D. Marshall, and B. O. Baumgaertner, Transient prophylaxis and multiple epidemic waves, AIMS Mathematics, 7 (2022), pp. 5616–5633, https://doi.org/10.3934/math.2022311.
  * [38] D. Weston, K. Hauck, and R. Amlôt, Infection prevention behaviour and infectious disease modelling: A review of the literature and recommendations for the future, BMC Public Health, 18 (2018), https://doi.org/10.1186/S12889-018-5223-1.