Physica D 448 (2023) 133714 Contents lists available at ScienceDirect Physica D journal homepage: www.elsevier.com/locate/physd Saddle–node bifurcation of limit cycles in an epidemic model with two levels of awareness ∗ David Juher, David Rojas, Joan Saldaña Departament d’Informàtica, Matemàtica Aplicada i Estadística, Universitat de Girona, Girona 17003, Catalonia, Spain article info a b s t r a c t Article history: In this paper we study the appearance of bifurcations of limit cycles in an epidemic model with two Received 4 October 2022 types of aware individuals. All the transition rates are constant except for the alerting decay rate of Received in revised form 20 January 2023 the most aware individuals and the rate of creation of the less aware individuals, which depend on the Accepted 20 February 2023 disease prevalence in a non-linear way. For the ODE model, the numerical computation of the limit Available online 11 March 2023 cycles and the study of their stability are made by means of the Poincaré map. Moreover, sufficient Communicated by M. Simpson conditions for the existence of an endemic equilibrium are also obtained. These conditions involve a Keywords: rather natural relationship between the transmissibility of the disease and that of awareness. Finally, Epidemic models stochastic simulations of the model under a very low rate of imported cases are used to confirm the Awareness scenarios of bistability (endemic equilibrium and limit cycle) observed in the solutions of the ODE Bifurcations model. Limit cycles © 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND Stochastic simulations license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction a decreasing function of the prevalence of the disease [3,5–7]. In this mean-field formulation of the incidence term, β depends on The role of human behaviour has been increasingly consid- the contact rate as well as on the probability of transmission dur- ered in epidemiological modelling since the early 2000s [1]. The ing an infectious contact. So, its reduction can reflect a diminution spread of COVID-19 has highlighted even more its important in the number of social contacts (social distancing), the adoption role in the progress of infectious diseases. Besides institutional of measures to prevent infection while keeping the same contact measures as mobility restrictions, mandatory use of facemasks, rate (decrease of the infection probability), or both. or school closings, self-initiated individual behaviours related On the other hand, it is well known that the perception of to risk aversion are recognized as a driving force in epidemic infection risk is uneven among susceptible individuals [8]. One dynamics [2,3]. Another important example of the impact of self- way to introduce some heterogeneity in risk-taking propensity initiated individual behaviours on epidemic spreading is given by has been to include more types of uninfected individuals char- acterized by their level of responsiveness to risk. For instance, the evolution of sexually transmitted diseases (STDs). The current the Susceptible–Aware–Infectious–Susceptible (SAIS) model con- resurgence in the number of cases of STD such as gonorrhea siders a new class of non-infected individuals with a higher and syphilis began in the mid-1990s, after the striking decline in risk aversion than the susceptible ones, the so-called aware or the number of STD cases following the appearance of the HIV in alerted individuals, who are characterized by a lower transmis- the early 1980s and the subsequent widespread use of condoms. sion rate [9]. However, a lower perception of risk following the introduction of A basic ingredient in such a modelling approach is the trans- antiretroviral therapies for HIV led to a decrease in condom use mission of awareness among individuals [10]. In [4] the authors and the current explosion of STD cases (see [4] and references considered an SAIS model where alerted individuals were able to therein). transmit awareness by convincing non-aware individuals to take One way to model such behavioural changes in deterministic preventive measures against the infection, which is an example models is to modify the incidence term β SI where β denotes of self-initiated individual behaviour. Moreover, a new class of the rate of disease transmission, and S and I are the number of aware individuals, the so-called unwilling (U) individuals, is also susceptible and infected individuals, respectively. The simplest introduced. They are characterized by a lower level of alertness way to modify it is by assuming that β is no longer constant but which is translated into a lack of willingness to transmit aware- ness to susceptible individuals. The existence of this second class ∗ Corresponding author. of aware individuals turns out to be necessary to have oscillatory E-mail addresses: david.juher@udg.edu (D. Juher), david.rojas@udg.edu solutions of the SAUIS model with no births and deaths in the (D. Rojas), joan.saldana@udg.edu (J. Saldaña). population. https://doi.org/10.1016/j.physd.2023.133714 0167-2789/© 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by- nc-nd/4.0/). D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 Precisely, we prove the existence of endemic equilibria for the model with ε = 0 and a general function σm (i) under biologi- cally meaningful conditions by proving and applying an extended version of the Poincaré–Miranda theorem on a triangular do- main. Moreover, for ε > 0, a set of parameters satisfying these conditions, and using m as a tuning parameter, we obtain two regions of bistability that appear from a saddle–node bifurcation of limit cycles, one for low values of m and the second one for high values of m. These two scenarios of bistability are confirmed by performing stochastic simulations of the transitions involved in the model. To our knowledge, this is the first time that such bifurcations are reported in epidemic models with dissemination of awareness, whether deterministic or stochastic. 2. SAUIS-ε model with varying coefficients Each individual in a population can be in one of the following four states: S (susceptible), A (aware), U (unwilling), and I (in- fected). The model assumes that aware individuals are created Fig. 1. Shape of function σm (i) with η = 0.4 for different values of m. at alerting rates αi and αa from susceptible ones after being in contact with infected and aware individuals, respectively. Aware individuals experience an alerting decay and become unwilling at a rate δa , while unwilling individuals also appear at rate νa The inflow of new susceptible individuals in the population is from contacts between susceptibles and aware individuals and a key factor in mean-field epidemic models to observe periodic they become susceptible at a rate δu . The infection transmission solutions [6]. In dynamic networks models, link dynamics can rates for susceptible, aware, and unwilling individuals are β , βa , also play this role [11]. However, even without demographic and βu , respectively, while the recovery rate from infection is δ . processes, behaviourally-induced epidemic oscillations can also Moreover, following the SAUIS-ε model introduced in [13], be expected to occur when individuals experience a decline in we consider the entry of imported cases (infections contracted awareness as a result of preventive measures taken over long abroad) at a very low rate ε > 0. This fact will prevent stochastic periods of time combined with low disease prevalence. The epidemic oscillations from extinction and, at the same time, the WHO report entitled ‘‘Pandemic Fatigue’’ defines this decline dynamical properties of the solutions remain close to those of as ‘‘demotivation to follow recommended protective behaviours, the deterministic ODE model with ε = 0. From a demographic emerging gradually over time and affected by a number of emo- point of view, it is assumed that the arrival of imported cases tions, experiences and perceptions’’ [12]. The existence of this is a sporadic event and, therefore, does not significantly alter type of oscillations was, indeed, proved in [4] by analysing the the prevalence of the disease. In stochastic simulations, it will occurrence of a Hopf bifurcation from an endemic equilibrium mean on average less than one imported case per mean infectious of the SAUIS model. Later, the robustness of such oscillations period when the population size is N = 105 because only was confirmed in [13] under the assumption of a low rate ε of non-infected individuals can bring infections from abroad. imported cases (infections contracted from abroad) by means of Finally, as explained in the introduction, we assume that the stochastic simulations on random networks. rate of awareness decay δa and the rate νa at which susceptible In this paper, we explore an extended version of the SAUIS-ε hosts become unwilling due to a contact with an aware host both model in [13] which considers that awareness dynamics changes depend on the fraction of infected individuals through a reduction abruptly when disease prevalence crosses a threshold value η. factor given by the function σm (i). Precisely, the rate νa of creation of unwilling individuals and Here it follows a summary of all transitions defining the the rate of awareness decay δa are modulated by the following SAUIS-ε model for a population of N individuals: function σm (i) of the fraction i of infected individuals: β0 βa0 βu0 I + S −→ 2I , I + A −→ 2I , I + U −→ 2I 1 σm (i) = , m ≥ 1. αi0 αa0 ν 0a (1) 1 + (i/η)m I + S −→ I + A, A + S −→ 2A, A + S −→ A + U δ δa δu ε The sharpness of the reduction of these two rates is controlled I −→ S , A −→ U , U −→ S , {S , A, U } −→ I by the parameter m, while η is the half-saturation constant (see Fig. 1). In particular, for m ≫ 1, 1 − σm (i) becomes closer to where δ a = δa σm (i) and the unit step function θ (i − η). Sudden behavioural changes can β βa βu αi occur, for instance, with the appearance of new cases of a highly β0 = , βa0 = , βu0 = , αi0 = , N −1 N −1 N −1 N −1 fatal infectious disease such as that caused by the Ebola virus, a (2) reemerging disease that generated panic among people during its αa νa σm (i) αa0 = , ν 0a = 2014 outbreak [14]. N −1 N −1 In contrast to other papers where awareness is considered in are the corresponding rates per single contact. terms of non-constant infection transmission rates (see, for in- After replacing N − 1 by N in (2) for simplicity, the resulting stance, [3,5–7,15]), here we will focus on the awareness dynamics ODE system governing the epidemic dynamics is given by: themselves and their role in the appearance of periodic solutions da (oscillatory epidemics). We are interested in how the behaviour = αi s i + αa s a − βa a i − δa σm (i) a − εa, βa < β, dt of solutions is affected by the reduction of both the decay of du awareness and the creation of unwilling individuals. Moreover, = δa σm (i) a + νa σm (i) s a − βu u i − δu u − ε u, βu < β, (3) dt we will use a parameter, m, which is not a transmission rate but di it is related to the behavioural response itself. = (β s + βa a + βu u − δ )i + (1 − i)ε, s + a + u + i = 1, dt 2 D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 δa (β−βu )−αa (δ−βu ) where s, a, u, and i denote the fractions of hosts in the S, A, U, which vanishes at a = 0 and a = αa (βu −βa ) . Elementary and I compartments, respectively. The differential equation for s β−δ computations show that this second root is greater than β−βa if has been omitted because it is redundant. and only if (β − βu )(δa (β − βa ) − αa (δ − βa )) > 0, 3. Equilibria which follows from the hypotheses. Thus f1 < 0 on the hy- The natural state space of system (3) is Ω := {(a, u, i) ∈ R : 3 potenuse (note that the coefficient of a2 is strictly positive since 0 ≤ a + u + i ≤ 1}. The existence of imported cases from abroad βu > βa ). β (β−β ) (β−δ )βu +βδu guarantees that the vector field defined by this system on the On the other hand, f2 (0, u) = u β u u2 − β u is β−δ boundary of Ω points strictly towards its interior. In particular, negative for all u ∈ (0, β−βu ) and vanishes at u = 0, and this implies that Ω is positively invariant under the flow defined by the solutions of system (3) and, moreover, the non-existence σm (i) ( f2 (a, 0) = −βa νa a2 + (βδa + δνa )a ) of disease-free equilibria for this model. β On the other hand, since σm (0) = 1, the same analysis of βδa +δνa vanishes at a = 0 and a = βa νa > 1, so f2 (a, 0) > 0 for all the bifurcations from the two disease-free equilibria (DFE) of β−δ the model with ε = 0 done in [13] works for our system. For a ∈ (0, β−βa ). instance, taking ε as a bifurcation parameter, it follows that one Therefore, f (a, u) := (f1 (a, u), f2 (a, u)) satisfies the assump- interior equilibrium of (3) comes from the bifurcation of a DFE tions of Theorem A.2 (see the Appendix) and there exists a of the system with ε = 0 when β < δ . Precisely, either the common zero of f1 (a, u) and f2 (a, u) inside the triangle formed DFE e∗1 = (0, 0, 0) enters Ω for ε > 0 if αa < δa , or the DFE by the plane (4) inside Ω , which corresponds to an endemic e∗2 = (a∗0 , u∗0 , 0) with a∗0 = δu (1 − δa /αa ) / (δa (1 + νa /αa ) + δu ) equilibrium. □ and u∗0 = δa /δu (1 + νa /αa ) a∗0 enters Ω for ε > 0 if αa > δa . αa β In both cases, the interior equilibrium of system (3) bifurcates Note that, if δa < δ and the first hypothesis of the lemma is β−x from an asymptotically stable DFE and is only maintained by fulfilled, then the second one is guaranteed because g(x) = δ−x the presence of imported cases. So, such an equilibrium is not a is an increasing function of x and, hence, we have proper endemic equilibrium because it does not result from the αa β β − βa disease transmission within the population. < < . For β > δ and taking βa as a bifurcation parameter, it δa δ δ − βa ∗ (follows that e2 is still ) asymptotically stable if βa < βac := β − In other words, if the basic reproduction number of the disease β − δ − (β − βu )u0 /a0 . In this case, an interior equilibrium fed ∗ ∗ in an awareness-free population, β/δ , is larger than 1 and, also, by the imported cases bifurcates from it. So, from now on we will is larger than that of the awareness transmission in a wholly assume that β > δ and βa > βac to guarantee that no interior susceptible population, αa /δa , then the existence of an endemic equilibrium for ε > 0 arises from a DFE with ε = 0 and, hence, equilibrium for ε = 0 is guaranteed. that any interior equilibrium corresponds to the perturbation of With respect to the endemic equilibria of system (3) for ε > 0, an endemic equilibrium of the system with ε = 0. we know that any equilibrium whose existence is guaranteed by Endemic equilibria are, in general, very difficult to determine Lemma 3.1 will persist inside Ω for ε > 0 small enough under analytically. When ε = 0 we can easily see that any endemic the classical transversal condition thanks to the implicit function equilibrium lies inside the plane theorem. δ βa βu ( ) ( ) 1− − 1− a− 1− u − i = 0. (4) 4. Saddle–node bifurcation of limit cycles β β β The following result gives sufficient conditions for the existence A bifurcation that passes from the exhibition of two hyper- of at least one endemic equilibrium point of the model (3) with bolic limit cycles of different stability in the phase portrait to ε = 0. The proof relies on a version of the Poincaré–Miranda the- the absence of such limit cycles through their collision in a orem in a triangular domain, which we include in the Appendix non-hyperbolic semistable limit cycle is called a saddle–node bi- for completeness. furcation of limit cycles. A classic scenario in the plane where this bifurcation appears is through a bistable configuration, when a Lemma 3.1. System (3) with ε = 0 has an endemic equilibrium stable equilibrium point is surrounded by an also stable limit β−β point in Ω if 0 ≤ βa < βu < δ < β and αδ a < δ−β a . cycle. Therefore a second limit cycle, in this case unstable, sepa- a a rates the basins of attraction of both stable objects. A continuous dependence on parameters may cause the two limit cycles col- Proof. Substituting (4) into the first and second equations of (3) lide and initial solutions that were attracted to the oscillatory we obtain two continuous functions in the variables (a, u), f1 motion are then attracted to the equilibrium. This bifurcation is and f2 , respectively. We find endemic equilibria in the common not exclusive of the plane, although for greater dimensions the zeros of f1 (a, u) and f2 (a, u). The intersection of the plane (4) with existence of the equilibrium is not required. Ω projected to the (a, u)-plane is the right triangle with vertex β−δ β−δ β−δ β−δ In order to find limit cycles and to determine their stability, (0, 0), ( β−β , 0) and (0, β−β ). Notice that β−β , β−β < 1 since a u a u the Poincaré map is used. In three-dimensional vector fields, as it βu , βa < δ . The hypotenuse of the triangle is given by substituting is the case for the model under consideration, the Poincaré map i = 0 in Eq. (4). (also known as first-return map) is a two-dimensional discrete α β (β−β ) α (ββ +βδ−2βu δ ) On the one hand, f1 (0, u) = i u β 2 u u2 − i u β 2 u+ map from a plane transversal to the flow of the system located αi δ (β−δ ) β−δ β−δ near the periodic orbit to itself. The image by the Poincaré map β2 is positive for u ∈ [0, β−βu ) and vanishes at u = β−βu . On of each point on the transversal plane is the next intersection the hypotenuse, f1 takes the value point of the flow on the plane. Fixed points of the Poincaré map β−δ−(β−βa )a a correspond to limit cycles and the stability of such fixed points f1 (a, β−βu ) = (αa (βu −βa )a +αa (δ −βu ) −δa (β −βu )) gives the stability of the periodic orbit. β − βu 3 D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 Fig. 2. Two solutions of system (3) for m = 3.5 tending to the limit cycle (in green) whose dynamics is close to the plane (4) (in light blue). 4.1. Numerical computation of limit cycles and their stability Table 1 Standing values of some parameters. β βa βu αi αa νa δ δa δu ε η Analytic treatment of bifurcations involving limit cycles are 2 0 0.5 0.001 0.015 3 1 0.01 0.03 10−5 0.4 only available when the expression of the periodic orbit is known as a function of the parameter producing the bifurcation. Since such expression is usually not computable in applications, nu- merical methods are the common technique for the detection of Then the Jacobian matrix of the Poincaré map can be written as saddle–node bifurcations of limit cycles. Dg(P(x))DϕT (x) (x) Consider a system of differential equations ẋ = f (x) for x ∈ Rn DP(x) = −f (P(x)) + DϕT (x) (x). Dg(P(x))f (P(x)) and the associated flow ϕt (x), t ∈ R. Consider a hypersurface Σ of Rn transversal to the vector field and assume there exists Finally, DF = DP − I and the iterative procedure (5) can be x0 ∈ Σ such that ϕT0 (x0 ) ∈ Σ for some minimal time T0 > used to locate limit cycles (we refer the reader to [16] for more 0. Continuous dependence on initial conditions of the system details). Once the limit cycle is located, the eigenvalues of the provides the existence of a neighbourhood U of x0 and a function monodromy matrix DϕT (x) (x) at the limit cycle, the so-called T : U → R such that ϕT (x) (x) ∈ Σ for all x ∈ U, the so-called Floquet characteristic multipliers, give the stability of the limit time-return map. The map P : Σ → Σ defined by P(x) = ϕT (x) (x) cycle found. is the Poincaré map or first-return map of the section Σ . Fixed points of P correspond to limit cycles of the system. So, 4.2. The saddle–node bifurcation of limit cycles in the model in order to locate them, the Newton–Raphson’s method can be applied to the distance function F (x) = P(x) − x near the limit The version of the SAUIS-ε model (3) exhibits a saddle–node cycle we wish to locate. The iterative procedure bifurcation of limit cycles in a certain region of the parameter space using m as a bifurcation parameter. Throughout the paper xn+1 = xn − DF (xn )−1 F (xn ), we take fixed values for some parameters. See Table 1. For those parameters Lemma 3.1 ensures the existence of one equilibrium, where DF denotes the Jacobian matrix of F and xn ∈ Σ , produces e∗ (m) = (a∗ (m), u∗ (m), i∗ (m)), which in this case is unique and successively better approximations of the periodic orbit. Rather lies inside the plane (4). For the computation of the Poincaré map than computing the inverse of the Jacobian matrix, the usual and we consider the plane Σ := {g(a, u, i) = a − a∗ (m) = 0} as a more numerically stable procedure is to solve the linear system Poincaré section near the limit cycles. For values of the parameter m < mh ≈ 3.698, the equilibrium DF (xn )(xn+1 − xn ) = −F (xn ). (5) e∗ (m) is unstable and the system presents an attractive limit Newton–Raphson’s method is an effective way to find isolated cycle that is mostly restricted to that plane (see Fig. 2). As the periodic orbits as long as the Jacobian matrix DF is non-singular. parameter increases, the stability of the equilibrium changes, pro- In order to integrate the solution we use Runge–Kutta’s method ducing a subcritical Hopf bifurcation at m = mh and a bistability RK45. At the same time the solution ϕt (x) is integrated, we also scenario for m > mh : the equilibrium and the limit cycle (see top integrate the variational equation panels in Fig. 3). As a consequence of the subcritical character of the Hopf bifurcation, an unstable limit cycle is born from the Ẏ = Df (ϕt (x))Y , Y (0) = In equilibrium. As long as the parameter continues increasing, both limit cycles (stable and unstable) start approaching each other. obtaining the monodromy matrix DϕT (x) (x). Therefore the differ- At the saddle–node bifurcation point m = msn1 ≈ 3.761, both ential of the Poincaré map can be computed as limit cycles collide in a non-hyperbolic semistable limit cycle and disappear, leaving the equilibrium as the only stable scenario DP(x) = f (P(x))DT (x) + DϕT (x) (x). (see bottom panels in Fig. 3). The dynamics remain similar for msn1 < m < msn2 and, at the second bifurcation point m = msn2 ≈ If the hypersurface Σ is defined by {g(x) = 0} and it is traversed 16.057, a second saddle–node bifurcation of limit cycles occurs. from {g(x) < 0} to {g(x) > 0}, by implicit derivation of g(x) = 0 In this case, a non-hyperbolic semistable limit cycle appears and we can find the differential of the time-return map splits into two hyperbolic (stable and unstable) limit cycles as Dg(P(x))DϕT (x) (x) m > msn2 , and this happens without any change of the stability DT (x) = − . of the equilibrium e∗ (m) which remains always asymptotically Dg(P(x))f (P(x)) 4 D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 Fig. 3. Projections on the (a, u)-plane of some solutions of system (3) for m = 3.69 (top left), m = 3.73 (top middle), m = 3.76 (top right), m = 3.76082 (bottom left), m = 3.77 (bottom middle) and m = 3.8 (bottom right). The orbit in green corresponds to the hyperbolic attractive limit cycle, the orbit in red to the hyperbolic repulsive limit cycle and the orbit in blue to the non-hyperbolic semistable limit cycle. Blue circles mark the initial conditions. Fig. 4. Projections on the (a, u)-plane of some solutions of system (3) for m = 16 (top left), m = 16.05 (top middle), m = 16.0575 (top right), m = 16.05755 (bottom left), m = 16.06 (bottom middle) and m = 17 (bottom right). The orbit in green corresponds to the hyperbolic attractive limit cycle, the orbit in red to the hyperbolic repulsive limit cycle and the orbit in blue to the non-hyperbolic semistable limit cycle. Blue circles mark the initial conditions. stable (see Fig. 4). The limit cycles separate each other until DϕT (x) (x) of the stable (solid blue line) and unstable (dashed red a position which is qualitatively unchanged as m increases. In line) orbits, showing the stability of each limit cycle. We point out Fig. 5 we show the previous described bifurcation phenomena, that, since DϕT (x) (x) is the monodromy matrix of a limit cycle, where the amplitudes of the stable limit cycle (solid line) and one of its eigenvalues is always 1 (the one with eigenvector unstable limit cycle (dashed line) with respect to the proportion orthogonal to the section Σ ). Moreover, since the motion is of infected individuals are displayed. The amplitude is computed rapidly almost captured by the plane (4), a second eigenvalue is as the difference between the largest and smallest value of the close to zero. The stability of the limit cycles is then given by the proportion of infected individuals along the orbit. Zero amplitude remaining third eigenvalue. On the left-hand panel we can see corresponds to the equilibrium e∗ (m). In Fig. 6 we represent how the unstable limit cycle appears for m = mh , with its third the Floquet characteristic multipliers of the monodromy matrix eigenvalue being larger than one. When approaching m = msn1 5 D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 Fig. 5. Bifurcation diagram for ε = 0 (left) and ε = 10−5 (right). Fig. 6. Floquet characteristic multipliers of the stable limit cycle (solid blue line) and the unstable limit cycle (dashed red line). the unstable eigenvalue tends to 1, as it does the third eigenvalue its weight. More precisely, we consider the interval [0, y] where of the stable limit cycle, producing the saddle–node bifurcation of limit cycles and the semi-stability of the orbit. On the right-hand y = IS β 0 + IAβa0 + IU βu0 + IS αi0 + AS αa0 panel we can see that at m = msn2 the semi-stable limit cycle appears giving birth to the stable and unstable limit cycles for + AS ν 0a + I δ + Aδ a + U δu + (N − I)ε. m > msn2 . We point out the strength of the instability in this case, as we also see on the bottom right panel of Fig. 4, where We think [0, y] as the total event window, partitioned in ten orbits are rapidly pushed away from the unstable limit cycle. subintervals corresponding to the ten possible types of reactions listed in (1), the length of each subinterval being a summand in the previous equality. Then, a number x in [0, y] is chosen 5. Stochastic simulations uniformly at random. The subinterval containing x tells us the type of single event that will occur. Once chosen, the event takes System (3) is a continuous model for the averaged propagation place and A, I, U and S = N − A − I − U are accordingly updated. of an SAUIS-ε epidemics among the individuals in a fully mixed Finally, we choose at random a number 0 < z < 1, increase the population. So, we need to perform stochastic simulations in continuous time by − log(z)/y (see [13,17]) and proceed to the order to assess in which sense and to which extent the invariant next iteration. objects exhibited by system (3) are found also in a discrete Given a population size N, a combination of model parame- context. ters, and an initial condition a(0), u(0), i(0), we run 50 indepen- dent simulations, each corresponding to a random distribution of 5.1. General simulation setup a(0)N, u(0)N, i(0)N and (1 − a(0) − u(0) − i(0))N individuals having respectively the initial states of aware, unwilling, infected and As usual in the setting of continuous-time stochastic simula- susceptible. For any experiment, we store the evolution of a(t), tions, we use the well-known Gillespie algorithm (GA in what u(t) and i(t) as three time series of equally-spaced points in the follows) [17], which was originally designed to simulate a fully interval [0, T ], where T is the maximum running continuous-time mixed chemically reacting system. We have a population of N of the simulation. Throughout the paper, the caption of each re- individuals and the algorithm keeps trace of the total numbers S, ported figure obtained by simulation includes the specification of A, I, U of susceptible, aware, infected and unwilling individuals. the values of N, a(0), u(0), i(0) and T . Recall that some parameters At each iteration, a single event is chosen at random according to have the standing values shown in Table 1. 6 D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 Fig. 7. On top, two outputs from GA with parameter m = 3.7 with N = 105 (dotted line) and N = 106 (solid line), and initial conditions a0 = 0.039359, i0 = 0.231104, u0 = 0.306049 (left), a0 = 0.039359, i0 = 0.331104, u0 = 0.172716 (centre), a0 = 0.039359, i0 = 0.431104, u0 = 0.039382 (right) and same initial random seed. On bottom, two outputs from GA with parameter m = 3.7 with N = 106 and same initial conditions as before with different initial random seed. 5.2. Detection of the bistability regime shown in Fig. 5. A more clear example is given in Fig. 8 for the parameter values m = 17 and m = 18. With the aim of producing an analogous stochastic version of The previous discussion motivates to choose a more useful the numerical bifurcation diagram depicted in Fig. 5 we follow representation of the different realizations (see [18] for a similar the procedure mentioned in the previous section. We perform 50 approach). To this end, and with the aim of producing a stochastic independent simulations for each value of the parameter m and bifurcation diagram, for each parameter m and each initial con- initial conditions starting from the location of the equilibrium dition we represent the maximum amplitude of each individual e∗ (m) = (a∗ (m), i∗ (m), u∗ (m)), keeping a(0) = a∗ (m) fixed, realization of the experiment in the time interval [1000, 10000]. increasing the value of i(0) from i∗ (m) by 0.01 and computing For instance, the maximum amplitude of approximately 1000 u(0) according to the equality (4). In this way, for each m ap- different realizations is shown in Fig. 9. We omit the first part of proximately 20 different initial conditions are considered inside the realization seeking for stationarity of the time series. On the a straight line lying on the plane (4). The larger i(0), the farther left panel, a transition can be noticed around the parameter value from the equilibrium the initial condition. m = 4, where amplitudes pass from 0.45 to nearly neglectable In stochastic epidemic models the average of the data is the (around 0.1). On the right panel, a similar situation occurs near usual way to construct a single signal to compare with the an- m = 16. These two values are close to the bifurcation points alytic dynamics. This is so because ODE systems are a good msn1 ≈ 3.761 and msn2 ≈ 16.057 numerically computed for approximation of the mean of realizations of stochastic processes the system (3). However, a clear difference on the abruptness of in systems with a large number of components. However, al- the transition can be appreciated between the two bifurcation though the standard mean of trajectories is a good option when points. On the right-hand panel, the bistability scenario is clearly dealing with high prevalence endemic equilibria, it does not represented and only a small amount of realizations near msn2 always work when a system exhibits fluctuating dynamics and differ from their amplitudes. However, in the left-hand panel even less in the presence of bistability. In the three lower panels the change in the amplitude seems to be more continuous. This of Fig. 7 we see two different realizations of the same experi- happens for two reasons. The first one is, as usual, related to the ment with identical initial conditions but different initial random number of individuals. In Fig. 10 a zoom near m = 3.8 is given, seed. Although they start identically, it is clear that stochasticity showing two panels with N = 105 (left) and N = 106 (right). reveals differences in the global behaviour beyond the expected We note that the region of bistability is narrower and closer to small perturbations. This phenomenon is due to the bistability of msn1 in the right-hand panel. Therefore we can expect a better the system, where small stochastic perturbations may change a diagram as N increases. The second reason is dynamical. First, the realization from one attraction basin to the other. Of course, the region of bistability in the parameter space for mh < m < msn1 greater the number of individuals N, the more similitude between is really small compared with m > msn2 . Second, the stability of realizations and with the deterministic scenario. The three top the limit cycles is also very different. For mh < m < msn1 the panels of Fig. 7 show how lower number of individuals produce unstable limit cycle has very weak repulsion as shown in the left more instability of the realization (for instance, abrupt changes of panel of Fig. 6. As a consequence, the dynamics near the unstable the amplitude) than higher number of individuals. The bistability limit cycle can be misunderstood as stochastic periodic solutions region in the previous figure is sensitive and very narrow, as due to the slow decay of the amplitude (even small stochastic 7 D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 Fig. 8. On top, two outputs from GA with parameter m = 17, N = 105 and initial conditions a0 = 0.033240, i0 = 0.437330, u0 = 0.039240 (left), a0 = 0.033240, i0 = 0.447330, u0 = 0.025900 (centre), a0 = 0.033240, i0 = 0.457330, u0 = 0.012570 (right). On bottom, a single output from GA with parameter m = 18, N = 105 and initial conditions a0 = 0.033240, i0 = 0.437420, u0 = 0.039140 (left), a0 = 0.033240, i0 = 0.447420, u0 = 0.025810 (centre), a0 = 0.033240, i0 = 0.457420, u0 = 0.012480 (right). Fig. 9. Maximum amplitude of each GA output in the time interval [1000, 10000] with N = 105 . perturbations may counter the decay). This is clearly represented addition of new compartments for alerted/responsive individuals in Fig. 3. On the other hand, for m > msn2 the repulsion of the that have reduced transmission rates [9–11,20,21]. unstable limit cycle is very strong as shown in the right-hand In this paper, we combine both approaches by considering an panel of Fig. 6 and this makes the realizations either stay on the epidemic model without demography which includes two types stable limit cycle or rapidly tend to the equilibrium, as shown in of aware individuals who are distinguished by their willingness Fig. 8. to convince susceptible individuals to adopt preventive measures. For the model with constant rates, we know that oscillatory 6. Conclusions solutions can appear as a consequence of a supercritical Hopf bifurcation from the endemic equilibrium [13]. Now, assuming The existence of a high number of infectious cases in a popula- that the rate of alerting decay δa as well as the rate of creation of tion during an epidemic can modify our individual behaviour and new unwilling individuals νa decrease by a nonlinear reduction how we relate to others. In turn, behavioural changes modify the factor σm (i) ∈ (0, 1] as the prevalence i of the disease increases, epidemic spread itself. This interplay has long been considered we have shown the existence of two scenarios where a bistable in many papers dealing with classic compartmental models with configuration with a stable limit cycle and a stable endemic incidence functions depending on the disease prevalence in a equilibrium occur. Precisely, we assume an abrupt change of both nonlinear way [6,7,15,19]. An alternative approach is based on the rates when disease prevalence crosses a threshold value η (the 8 D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 Fig. 10. Maximum amplitude of each GA output in the time interval [1000, 10000] with N = 105 (left) and N = 106 (right). half-saturation constant). Below this threshold, the value of the nonlinear rates that depend on the prevalence level in a popula- rates are close to their maximum values (δa and νa ), whereas tion. This has been obtained under a choice of parameters values they clearly decrease above it. The sharpness of this change is which assumes a much faster transmission among susceptible controlled by a parameter m which determines the slope of σm (i) individuals of both infections and low level of awareness than at i = η (σm′ (η) = −m/(4η)). In both scenarios, the parameters the creation of fully aware individuals. These dynamics are also values are in agreement with the sufficient conditions we have observed in stochastic simulations of the model assuming a (very obtained for the existence of, at least, one endemic equilibrium. low) rate of imported cases with large enough population sizes. Values of m ≫ 1 can be associated with radical changes in Note that the model does not include vaccination. Therefore, the self-initiated individual behaviour when the prevalence level behavioural aspects such as vaccine hesitancy, a delay in accept- is close to η. In this case, σm (i) ≈ 0 for η < i ≤ 1 which ing or refusing vaccines despite the availability of vaccination implies almost no decay of awareness and almost no creation of services [23], have not been considered in this paper. unwilling individuals (only fully aware individuals are created). For the parameters considered in the paper, a bistable configu- CRediT authorship contribution statement ration is always the case for m > 16.06 (σm′ (η) < −10) after the occurrence of a saddle–node bifurcation of limit cycles. This David Juher: Conceptualization, Methodology, Software, Writ- configuration is clearly observed in the stochastic simulations of ing – original draft. David Rojas: Conceptualization, Formal anal- the epidemic process with a very low rate of imported cases due ysis, Methodology, Software, Writing – original draft. Joan Sal- to the strongly repulsive character of the unstable limit cycle daña: Conceptualization, Supervision, Methodology, Writing – lying between the stable one and the endemic equilibrium. For original draft. lower values of m (here m < 4), the reduction of both rates is not so abrupt and bistability is only present for a narrow range Declaration of competing interest of values (3.698 < m < 3.761 with σm′ (η) ≈ −2.3) once a subcritical Hopf bifurcation has occurred. In this case, σm (i) is The authors declare that they have no known competing finan- clearly positive for η < i ≤ 1. For m < 3.698, the smoothness cial interests or personal relationships that could have appeared of the transition between high and low values of σm (σm′ (η) ∈ to influence the work reported in this paper. [−2.3, −0.625]) as well as the lower reduction of the two rates for i ≈ 1 make the endemic equilibrium unstable and allow for a Data availability stable limit cycle. The existence of imported cases (at a rate ε ) assumed in No data was used for the research described in the article. the present work has also been considered elsewhere (see, for instance, [18,22]). In addition to its suitability when modelling Acknowledgements epidemics in non-isolated populations, it prevents the stochas- tic extinction of oscillatory epidemics when disease prevalence This work is supported by the grants PID2019-104437GB-I00 reaches very low levels. From a deterministic point of view, as and PID2020-118281GB-C31 funded by MCIN/AEI/, Spain long as ε is small enough, the continuous dependence of solutions 10.13039/501100011033. The authors are members of the Con- on parameter values guarantees that the attractors of the model solidated Research Groups 2021 SGR 00113 of the Generalitat de with and without imported cases will be very close to each other Catalunya. D.R. is a Serra Húnter Fellow. (see Fig. 5). Higher values of ε are not expected to affect the main dynamical features of the model as long as the entry of Appendix imported cases is occasional (see, for instance, [13] where a value of ε one order of magnitude higher is also considered and it does In this Appendix we present a modified version of the not significantly alter the oscillatory dynamics of the model with Poincaré–Miranda theorem on the plane, which is an extended constant rates). version of the classical Bolzano’s theorem in higher dimension. In summary, we have shown the existence of bifurcations of Up to the authors’ knowledge, Poincaré–Miranda theorem is not limit cycles in the SAUIS model (without demography) when epi- trivially deduced in a planar triangular domain. We credit and demic spread and awareness transmission are coupled through thank professor Rafael Ortega for the idea of the proof, which relies on the following result of degree theory that we include 9 D. Juher, D. Rojas and J. Saldaña Physica D 448 (2023) 133714 Fig. 11. On the left, structure and assumptions of the Poincaré–Miranda theorem for triangles in Theorem A.2. On the right, a sketch of the behaviour of the function f (Γ (t)) for t ∈ [0, 1]. for the sake of completeness (see the Appendix on degree theory [3] J.S. Weitz, S.W. Park, C. Eksin, J. Dushoff, Awareness-driven behavior in [24].) changes can shift the shape of epidemics away from peaks and toward plateaus, shoulders, and oscillations, Proc. Natl. Acad. Sci. USA 117 (2020) Let Γ be a Jordan curve in R2 and let Ω be the open set 32764–32771. enclosed by Γ . Let f : Ω̄ → R2 be continuous such that f (x) ̸ = 0 [4] W. Just, J. Saldaña, Y. Xin, Oscillations in epidemic models with spread of for all x ∈ ∂ Ω = Γ . The degree of f in Ω , deg(f , Ω ), can be awareness, J. Math. Biol. 76 (2018) 1027–1057. computed as the winding number of f (Γ ) around the origin. [5] V. Capasso, G. Serio, A generalization of the Kermack-McKendrick deterministic epidemic model, Math. Biosc. 42 (1978) 43–61. [6] S. Ruan, W. Wang, Dynamical behavior of an epidemic model with a Theorem A.1. Let f : Ω̄ → R2 be continuous such that f (x) ̸ = 0 nonlinear incidence rate, J. Differential Equations 188 (2003) 135–163. for all x ∈ ∂ Ω . If deg(f , Ω ) ̸ = 0 then f has a zero in Ω . [7] M.E. Alexander, S.M. Moghadas, Periodicity in an epidemic model with a generalized non-linear incidence, Math. Biosc. 189 (2004) 75–96. Now we state the version of Poincaré–Miranda theorem on a [8] B. Guenther, M.M. Galizzi, J.G. Sanders, Heterogeneity in risk-taking during triangular domain, which can be easily generalized for any Jordan the COVID-19 pandemic: Evidence from the UK lockdown, Front. Psychol. 12 (2021) 643653. curve with similar assumptions. [9] F.D. Sahneh, F.N. Chowdhury, C.M. Scoglio, On the existence of a threshold for preventive behavioral responses to suppress epidemic spreading, Sci. Theorem A.2. Let Ω be a triangle and let f : R2 → R2 be Rep. 2 (2012) 632, http://dx.doi.org/10.1038/srep00632. a continuous function, f (x, y) = (f1 (x, y), f2 (x, y)). Consider the [10] S. Funk, E. Gilad, C. Watkins, V.A.A. Jansen, The spread of awareness and its impact on epidemic outbreaks, Proc. Natl. Acad. Sci. USA 21 (2009) boundary of Ω positively oriented and three distinguished points α , 6872–6877. β and γ as shown in Fig. 11. If f2 > 0 from α to β , f1 < 0 from β [11] A. Szabó-Solticzky, L. Berthouze, I.Z. Kiss, P.L. Simon, Oscillating epidemics to γ and f1 > 0 and f2 < 0 from γ to α , then f (x, y) has at least in a dynamic network model: stochastic and mean-field analysis, J. Math. one zero in Ω . Biol. 72 (2016) 1153–1176. [12] Pandemic Fatigue: Reinvigorating the Public To Prevent COVID-19, Docu- ment number: WHO/EURO:2020-1160-40906-55390, World Health Orga- Proof. Let us assume, with the aim of reaching contradiction, nization, 2020, https://apps.who.int/iris/bitstream/handle/10665/335820/ that f does not vanish on Ω . Let Γ : [0, 1] → R2 be a curve WHO-EURO-2020-1160-40906-55390-eng.pdf. travelling the boundary of Ω with Γ (0) = Γ (1) = α . Notice that, [13] D. Juher, D. Rojas, J. Saldaña, Robustness of behaviourally-induced oscilla- tions in epidemic models under a low rate of imported cases, Phys. Rev. by hypothesis, f does not vanish on the boundary of Ω . So, we E 102 (2020) 052301. can consider a continuous argument function θ : [0, 1] → R such [14] K.K.W. To, J.F.W. Chan, A.K.L. Tsang, V.C.C. Cheng, K.-Y. Yuen, Ebola virus that disease: a highly fatal infectious disease reemerging in West Africa, Microb. Infect. 17 (2015) 84–97. f (Γ (t)) = |f (Γ (t))| (cos θ (t), sin θ (t)) , [15] W. Liu, H.W. Hethcote, S.A. Levin, Dynamical behavior of epidemiological models with nonlinear incidence rates, J. Math. Biol. 25 (1987) 359–380. which is unique up to an additive constant 2π . By the hypotheses [16] J.M. Mondelo, Computing invariant manifolds for libration point missions, of the theorem, f (Γ (t)) lies on the first and second quadrants in: G. Baù, et al. (Eds.), Satellite Dynamics and Space Missions, in: Springer when travelling from α to β , and f (β ) lies on the second quadrant. INdAM Series, 34, pp. 159–223. [17] D.T. Gillespie, Stochastic simulation of chemical kinetics, Ann. Rev. Phys. Similarly, from β to γ , f (Γ (t)) lies on the second and third Chem. 58 (2007) 35–55. quadrants, finishing on the negative ordinate semi-axis. Finally, [18] M. Aguiar, B. Kooi, J. Martins, N. Stollenwerk, Scaling of stochasticity in from γ to α , f (Γ (t)) lies on the fourth quadrant ending on the dengue hemorrhagic fever epidemics, Math. Model. Nat. Phenom. 7 (2012) 1–11. positive abscissa semi-axis where it started (see Fig. 11). There- θ (1)−θ (0) [19] W.M. Liu, S.A. Levin, Y. Iwasa, Influence of nonlinear incidence rates upon fore, the winding number of the curve f (Γ (t)) is 2π = 1. the behavior of SIRS epidemiological models, J. Math. Biol. 23 (1986) By the definition of the degree, deg(f , Ω , 0) = 1 and then, by 187–204. Theorem A.1, f has at least one zero in Ω . □ [20] S. Funk, E. Gilad, V.A.A. Jansen, Endemic disease, awareness, and local behavioral response, J. Theoret. Biol. 264 (2010) 501–509. [21] I.Z. Kiss, J. Cassell, M. Recker, P.L. Simon, The impact of information References transmission on epidemic outbreaks, Math. Biosci. 225 (2010) 1–10. [22] R.M. Almarashi, C.C. McCluskey, The effect of immigration of infectives on [1] S. Funk, M. Salath, V.A.A. Jansen, Modelling the influence of human disease-free equilibria, J. Math. Biol. 79 (2019) 1015–1028. behaviour on the spread of infection diseases: a review, J. R. Soc. Interface [23] N.E. MacDonald, SAGE working group on vaccine hesitancy. Vac- 7 (2010) 1247–1256. cine hesitancy: Definition, scope and determinants, Vaccine 33 (2015) [2] S. Manrubia, D.H. Zanette, Individual risk-aversion responses tune epi- 4161–4164. demics to critical transmissibility (R = 1), R. Soc. Open Sci. 9 (2022) [24] R. Ortega, Periodic Differential Equations in the Plane: A Topological 211667. Perspective, De Gruyter, Berlin, Boston, 2019. 10