This is a preprint of: “Robustness of behaviourally-induced oscillations in epidemic models under a low rate of imported cases”, David Juher, David Rojas, Joan Saldaña, Phys. Rev. E, vol. 102(5), 052301: 1–15, 2020. DOI: [10.1103/PhysRevE.102.052301] Robustness of behaviourally-induced oscillations in epidemic models under a low rate of imported cases David Juher,∗ David Rojas,† and Joan Saldaña‡ Departament d’Informàtica, Matemàtica Aplicada i Estadı́stica, Universitat de Girona, Catalonia, Spain This paper is concerned with the robustness of the sustained oscillations predicted by an epidemic ODE model defined on contact networks. The model incorporates the spread of awareness among individuals and, moreover, a small inflow of imported cases. These cases prevent stochastic extinc- tions when we simulate the epidemics and, hence, they allow to check whether the average dynamics for the fraction of infected individuals are accurately predicted by the ODE model. Stochastic sim- ulations confirm the existence of sustained oscillations for different types of random networks, with a sharp transition from a non-oscillatory asymptotic regime to a periodic one as the alerting rate of susceptible individuals increases from very small values. This abrupt transition to periodic epi- demics of high amplitude is quite accurately predicted by the Hopf-bifurcation curve computed from the ODE model using the alerting rate and the infection transmission rate for aware individuals as tuning parameters. Keywords: epidemic models, awareness, oscillations, stochastic simulations. I. INTRODUCTION compartments. Other models assume one or more addi- tional compartments consisting of aware (A) individuals The importance of the interplay between epidemic [12–15]. Some of these works consider several awareness spreading and preventive behavioural responses in a glob- levels resulting from the assumption of a degradation in alized world has long been recognized and was specially the quality of the information as it is passed from one highlighted after the SARS outbreak of 2003 [1, 2]. The individual to another [16, 17]. In all these examples, the rise of the incidence rate of sexually transmitted diseases effect of preventive behaviour is to modify the values of (STDs) [3] and the current resurgence of measles [4] are the epidemic parameters like the probability of infection also examples of such an interplay. For STDs, increas- or the recovery rate. ing high risk sexual behaviour and novel sexual networks When the contact network structure of a population is are among factors responsible for their re-emergence, explicitly considered, the effect of behavioural responses whereas vaccine hesitance and distrust in public health can also affect the contact structure itself when mod- intervention programs are among behavioural factors re- elling social avoidance behaviours. This reduction of the sponsible for the rise of diseases like measles. exposure to disease has been modelled by means of pre- Risk perception is an important determinant of self- ventive disconnection from infectious neighbours [18–25] initiated, voluntary protective behaviour [5]. It consti- or, also, by replacing some infected nodes by healthy ones tutes the basic ingredient in many epidemic models to en- [26] leading, in both cases, to dynamical networks. Here capsulate human behaviour in their formulation [2]. For the assessment of disease prevalence is based on the in- instance, some extensions of classic deterministic com- dividual neighbourhood (local contacts), in contrast to partmental models include the impact of behaviour on homogeneous compartmental models where information disease transmission by assuming more general incidence about the prevalence is assumed to be globally available rates than the standard one (bilinear). The latter is pro- [2]. portional to the product of the number of susceptible (S) and infectious individuals (I), βSI, whereas its general- Under the previous modelling approaches, protective izations assume a saturation with respect to the number behavioural responses are triggered by the disease preva- of infectives in order to model a reduction of the contact lence. As long as these responses are based on the global rate in the presence of a high disease prevalence [6–9]. prevalence, one expects the likelihood of epidemic oscilla- Other model extensions take into account awareness tions to be high. Moreover, the linear dependence of the transmission among individuals. For example, one of standard incidence rates on the number of infected indi- them [10, 11] divides each epidemic compartment (S, I viduals implies that, when these oscillatory solutions oc- and R (recovered/immune)) into two subcompartments cur, they should pass through low prevalence levels due to of aware and unaware individuals, respectively, and in- the lack of an abrupt switching behaviour. A low preva- troduces the corresponding transition rates between sub- lence, in turn, will drive the number of aware individuals down as a consequence of a lower perception of the con- tagion risk, and the cycle repeats again with a new rise in the number of infectives. In fact, several ODE epidemic ∗ david.juher@udg.edu models with transmission of awareness [17, 27] or assum- † david.rojas@udg.edu ing self-initiated, voluntary vaccination [28] exhibit such ‡ joan.saldana@udg.edu periodic solutions under some values of the parameters. 2 In particular, the variability in the propensity of aware 2014-2016. However, two years later, this status was lost individuals to further propagate awareness as the driving after 991 confirmed cases in England and Wales in 2018, mechanism for sustained oscillations was proved in [17]. more than three times the number of cases (284) in 2017 More precisely, it was proved there that, in the absence [30]. Other European countries that have also lost the of demographics, standard incidence rates of infections WHO eradication status for measles in 2019 are Alba- and awareness do not lead to oscillations even though nia, Greece, and the Czech Republic [31]. Of course, the aware individuals disseminate awareness among suscep- role of imported cases is also important in new emergent tible ones, unless a class of individuals with a lower level diseases like Covid-19 where vaccination coverage is not of awareness (labelled U for unwilling to disseminate) is present at all. An illustrative example of their impor- also assumed. tance is given by the data of Covid-19 in Iceland (as of In this paper we study the robustness of the determin- June 3, 2020): 343 cases out of 1806 were reported to be istic oscillations of a model (which we call SAUIS-ε) that infected abroad [32]. is an extension of the SAUIS model studied in [17]. The Here, we consider the occasional introduction of new reason for this choice is twofold. First, because the only cases into the population, at a very small rate ε, the driving mechanism for the existence of periodic solutions rate of imported cases. From a mathematical point of in the SAUIS model is the variability in the propensity of view, these new cases will not destroy the deterministic alerted individuals to further propagate awareness in the oscillations obtained in [17] from a Hopf-bifurcation as population. Second, because the model does not con- long as ε is low enough and, at the same time, they pre- sider recruitment of susceptible individuals in terms of vent stochastic epidemic oscillations from extinction, as newborns or in terms of a prevalence-dependent recruit- it happens when periodic solutions reach very low levels ment of them into a core group, which turns out to be of prevalence (see Subsection IV C). an essential requirement for having sustained oscillations The main goal is, then, to analyse the SAUIS-ε model in many epidemic models with and without vaccination and to perform stochastic simulations of epidemics on [9, 28, 29]. networks in order to compare the regions of the parame- Another behavioural mechanism also responsible for ter space where oscillations occur under both (determin- the occurrence of periodic solutions in epidemic models istic and stochastic) approaches. Such a comparison will is preventive rewiring in adaptive contact networks. Such show that a significant part of the original region of the a mechanism aims to minimise the infection risk of sus- parameter space where deterministic oscillations occur is ceptible individuals while maintaining them connected to preserved when stochastic simulations are performed un- the network [18, 20, 21, 27]. So, both approaches try to der the presence of a low rate of imported cases. More- explain epidemic oscillations on a pure behavioural basis. over, we will see that, as expected, a low level of aware- On the other hand, damped epidemic oscillations have ness in the population at the moment of their introduc- been predicted in [16], where an SIR epidemic model tion raises the chances of sparking sustained transmission without demography is coupled with a sophisticated and, hence, of generating new waves of infection. This mechanism of degradation of the information quality. could be, for instance, the situation in many countries Such degradation is translated into individuals with dif- once lockdown restrictions for Covid-19 have been lifted ferent levels of awareness in the population. However, because most of their inhabitants are not immunized yet, the eventual depletion of susceptible nodes prevents the as it happened in Singapore with a second wave of Covid- occurrence of sustained oscillations. 19 infections within its population of foreigner workers In countries where an infectious disease has been de- [33]. clared eradicated, new cases can still occur, but these will be isolated and will have limited spread within the community as long as vaccination coverage is high II. THE SAUIS-ε MODEL ON A GENERAL enough. These new cases are usually imported either NETWORK from tourists, foreign workers, etc., or by local individu- als that have been infected abroad in regions where the According to the SAUIS model introduced in [17], each disease is endemic or large outbreaks are taking place. node of a network of size N can be in one of the follow- However, if the vaccination coverage in these countries ing four states: S (susceptible), A (aware), U (unwill- decreases because of lower levels of awareness (for in- ing), and I (infected). Given any node n, let us denote stance, in countries where vaccination is not mandatory), the probabilities that n will have the respective states such few cases can mount into major outbreaks. at time t by S n (t), An (t), U n (t) and I n (t). Note that In fact, regions that have achieved the WHO eradica- S n (t) + An (t) + U n (t) + I n (t) = 1 for all time t ≥ 0, since tion status for a given disease can lose it as a consequence each node has to be in one of the four states. Moreover, of a marked increase in the number of confirmed cases. for the sake of simplicity, we assume along the paper that For measles, a country attains this status when there the alerting rates per link αa0 , νa0 , αi0 , the infection trans- is no endemic transmission for 12 months in a specific mission rates per link β 0 , βa0 , βu0 , the decay rates δa , δu , geographic area. For instance, the UK achieved WHO and the recovery rate δ are the same for all the nodes. measles eradication status in 2017, based on data from For instance, if at time t a node n is aware and one of its 3 neighbours, m, is susceptible, then the probability that n The original SAUIS model in [17] contemplates the successfully alerts m during the time interval (t, t + ∆t) possibility that an infected individual, after recovering, is αa ∆t + o(∆t) provided that αa ∆t < 1. Similarly, the may become aware with probability p, unwilling with probability that a non-infected node contracts the infec- probability q, as well as susceptible with probability tion from abroad (imported case) during a time interval 1 − p − q. Under these two additional transitions, peri- of length ∆t is ε∆t+o(∆t). Here it follows a summary of odic solutions are also possible (cf. Figs. 7 and 8 in [17]) all transitions (or reactions) defining the SAUIS-ε model: but, since our ultimate goal is to provide evidence of ro- bustness of the oscillatory regime in non-deterministic epidemics and to simplify the analysis, we will assume β0 β0 β0 a I + S −→ I + I, I + A −→ u I + I, I + U −→ I +I p = q = 0 along the paper. For the sake of simplicity, we α0 i a α0 a ν0 have omitted them in the previous description of possible I + S −→ I + A, A + S −→ A + A, A + S −→ A+U transitions. δ a δ u δ ε I −→ S, A −→ U, U −→ S, {S, A, U } −→ I Let us derive the approximate discrete-time equations for the evolution of these probabilities. An exact (but un- feasible) description would require the probability of the Let us recall that the SAUIS model tries to account system being in any of the 4N possible states. So, to de- for the degradation of information quality among indi- rive approximate equations we will assume, as usual, that viduals that are aware of the epidemic situation. As the joint probability for nodes n and m to be respectively well as in the standard SAIS models, the transition in states X and Y is independent of the neighbourhood’s I + S → I + A represents the creation of a new aware in- configuration of n and m. That is, it equals the product dividual that has acquired first-hand information about of both probabilities. This hypothesis will allow us to the epidemic by means of a direct contact. Also, the tran- close the system without considering higher order terms sition A + S → A + A creates new aware individuals that for the joint probabilities (see [13] for a related discussion get indirect information from their acquaintances. Such for the SAIS model). new aware individuals have the same responsiveness as Let ∆t > 0 be small enough in such a way that, for the information disseminators, which is not always the every occurrence rate κ of a single event, the probability case. That is why the SAUIS model also includes the for this event to happen in the time interval (t, t + ∆t) A + S → A + U transition, where U stands for unwilling is κ∆t + o(∆t). For 1 ≤ n, m ≤ N , let anm be the to disseminate information. So, an unwilling individual (n, m) element of the N × N adjacency matrix of the has a lower level of awareness, in the sense that he or contact network, i.e. anm = 1, if the nodes n and m she does not try to convince other people about the risk, are first neighbours, and anm = 0 otherwise. With these having in addition a weaker behavioural response. ingredients, we can now write Q  An (t + dt) = An (t) 1 − δa dt − εdt − (1 − m (1 − anm βa0 dtI m (t))) Q Q  + S n (t) 1 − m (1 − anm αi0 dtI m (t)) m (1 − anm αa0 dtAm (t)) . The term multiplying An (t) corresponds to the event plicity, the probability of this third event is computed as that the node n keeps being aware at time t+dt provided 1 minus the probability that none of such neighbours suc- it was aware at time t, so it is 1 minus the sum of the ceeds. Analogously, the term multiplying S n (t) accounts probabilities of the three competing events that change for the probability that n is aware at time t + dt pro- the state A to another one: A → U with probability δa dt, vided it was susceptible at time t. Now observe that, ne- A → I with probability εdt, and the event that one or glecting Q terms of order o(dt), the expressions P ofnm the form nm m several infected neighbours of n succeed in performing the m (1 − a κdtX (t)) read as 1 − κdt m a X m (t). transition I +A → I +I (with probability βa0 dt). For sim- This yields P An (t + dt) = An (t)(1 − δa dt − εdt − βa0 dt m anm I m (t)) P P + S n (t)(1 − (1 − αi0 dt m anm I m (t))(1 − αa0 dt m anm Am (t))) P P P = An (t)(1 − δa dt − εdt − βa0 dt m anm I m (t)) + S n (t)(αi0 dt m anm I m (t) + αa0 dt m anm Am (t)), where in the second equality we have neglected again the terms of order o(dt). Subtracting An (t) to both sides of 4 the previous equation, dividing them by dt, and letting the time evolution for An (t). Proceeding along the same dt → 0, we obtain the differential equation governing lines for the other probabilities, we finally arrive at the following system of 3N ODEs: XN XN dAn (t) = anm (αa0 Am (t) + αi0 I m (t))S n (t) − βa0 anm An (t)I m (t) − (δa + ε)An (t), dt m=1 m=1 XN XN d U n (t) = δa An (t) + νa0 anm S n (t)Am (t) − βu0 anm U n (t)I m (t) − (δu + ε)U n (t), (1) dt m=1 m=1 XN dI n (t)  = anm β 0 S n (t) + βa0 An (t) + βu0 U n (t) I m (t) − δI n (t) + ε (1 − I n (t)) , dt m=1 P nm m where the equation for S n (t) is omitted because it is re- and ma I (t) = k i(t) because each node has the dundant (the sum of the nodal probabilities is always same degree k. So, the time evolution of these probabil- equal to 1). ities satisfies the following initial value problem (IVP): From the solution of system (1) endowed with an ini- tial condition, we can compute the expected number of da = αi s i + αa s a − βa a i − δa a − εa, βa < β, aware, unwilling and infectious nodes at time t by sum- dt ming the corresponding probabilities P n over the whole P net- du work, that is, N (t) = A (t), N (t) = U n (t), = δa a + νa s a − βu u i − δu u − εu, βu < β, (2) P A n U n dt NI (t) = n I n (t), and NS (t) = N − NA (t) − NU (t) − di NI (t). = (β s + βa a + βu u − δ)i + (1 − i)ε, dt Similar approaches to derive a system of equations for the probabilities for a node of being in one of several s + a + u + i = 1, endowed with the initial condition disease states have been previously introduced for the a(0) = a0 , u(0) = u0 , and i(0) = i0 . Here all the alerting study of epidemics on networks and have received dif- and infection transmission rates are per node and not per ferent names like, for instance, Microscopic Monte Carlo link, that is, αi = k αi0 , αa = k αa0 , νa = k νa0 , β = kβ 0 , Approach in a discrete-time setting [34], or N-intertwined βa = kβa0 , and βu = kβu0 . model in a continuous-time setting [35]. An extension of Note that an initial condition a(0) = a0 , u(0) = u0 , the latter to multilayer networks is given in [36]. i(0) = i0 of system (2) corresponds to the uniform ini- tial condition An (0) = a0 , U n (0) = u0 , I n (0) = i0 for 1 ≤ n ≤ N of system (1). Then, the local existence and III. THE SAUIS-ε MODEL ON REGULAR uniqueness of solutions for both IVPs implies that a so- RANDOM NETWORKS lution for (2) is also a solution (uniform by construction) for (1). That is, they are equivalent formulations of the To analyse system (1) we start by the simplest case. epidemic on regular networks (see Lemma 3.1 in [13] for So, let us consider the model over a random regular net- a proof of a similar result for the SAIS model). work (not necessarily fully connected) of degree k. As we On the other hand, notice that NA (t) = a(t)N , will see, in this particular case the solutions of the system NU (t) = u(t)N , and NI (t) = i(t)N . So, instead of think- (1) of 3N ODEs can be identified with the solutions of a ing of nodal probabilities, we can consider the expected much simpler system of three ODEs. fraction of aware, unwilling and infected individuals as On this sort of networks, every node has the same vul- the macroscopic description (state variables) of our sys- nerability against the disease (in this setting, the degree tem which is more convenient to compare theoretical pre- is the only characteristic that distinguishes one node from dictions with the outputs of the stochastic simulations of another). So, it is reasonable to assume the same ini- the epidemics. tial probabilities of being aware, An (0) = a0 , unwilling, As an incidental remark, it turns out that uniform U n (0) = u0 , and infected, I n (0) = i0 , for all nodes (what equilibrium solutions for the complete system (1) of 3N we call a uniform initial condition). equations are only possible over regular networks. To see Now we focus on uniform solutions of system (1), it, set An = A, U n = U , I n = I for all 1 ≤ n ≤ N defined as those solutions An (t), U n (t), I n (t) that are and denote the degree of any node n by kn . At the independent from n. So, we can write An (t) = a(t), equilibrium, we have that 0 = dAn (t)/dt = kn (αa0 A + U n (t) = u(t) and I n (t) = i(t) forP all 1 ≤ n ≤ N . Then, αi0 I)S − βa0 kn AI − (δa + ε)A. Considering a pair of differ- the sums in system (1) reduce to m anm Am (t) = k a(t) ent nodes n, m, from dAn (t)/dt = dAm (t)/dt we easily 5 get that kn = km . In consequence, all nodes must have and endemic equilibria P3 = (a∗ , u∗ , i∗ ) ∈ Ω with the same degree.     ∗ βa ∗ βu δ i =1− 1− a − 1− u∗ − > 0. β β β A. Behaviour of the equilibria Note that several distinct endemic equilibria may coex- Similarly to the SAUIS model in [17], the tetrahedron ist. When ε > 0 in the SAUIS-ε model, the first two Ω := {(a, u, i) ∈ R3 : 0 ≤ a + u + i ≤ 1} is positively kinds of equilibria cannot exist and all possible equilibria invariant under the flow of system (2). In fact, the vector are endemic. In fact, if ε is a small parameter the equi- field on the boundary of Ω points strictly towards its libria of SAUIS-ε can be interpreted as perturbations of interior for ε > 0. In consequence, there are no equilibria the equilibria P1 , P2 and P3 of the unperturbed SAUIS on the boundary of Ω. model. To study how these equilibria behave when ε When ε = 0 the SAUIS model may have three different is included in the model, we denote by fj (a, u, i; ε) the kinds of equilibria in Ω: the trivial equilibrium P1 = right-hand side of the jth equation in (2). Any equilib- (0, 0, 0), the disease-free equilibrium P2 = (a∗0 , u∗0 , 0) with rium e∗ (ε) = (a∗ (ε), u∗ (ε), i∗ (ε)) of the system is implic-      itly given by the equations fj (e∗ (ε); ε) = 0, j = 1, 2, 3. δu 1 − αδaa δa 1 − αδaa 1 + ανaa We can derive implicitly the equations with respect to ε a∗0 =   , u∗0 =   , and get νa δa 1 + αa + δu δa 1 + ανaa + δu ∂fj ∗ da∗ ∂fj ∗ du∗ ∂fj ∗ di∗ ∂fj ∗ (e (ε); ε) + (e (ε); ε) + (e (ε); ε) + (e (ε); ε) = 0, j = 1, 2, 3. ∂a dε ∂u dε ∂i dε ∂ε In the particular case of the equilibrium e∗ (0) = P1 = the introduction of new infection cases at a rate ε rather (0, 0, 0) the previous system for ε = 0 can be solved and than by the epidemic propagation itself. When the perturbation is considered from the equi- da∗ αi librium e∗ (0) = P2 = (a∗0 , u∗0 , 0), the previous equations = , dε ε=0 (β − δ)(αa − δa ) imply du∗ αi (δa + νa ) di∗ −1 = , = . dε ε=0 δu (β − δ)(αa − δa ) dε β − δ − (β − βa )a∗0 − (β − βu )u∗0 ε=0 di∗ 1 =− . Thus the condition such that the perturbation of the dε ε=0 β−δ disease-free equilibria is inside the region Ω for ε > 0 In order that the trivial equilibrium P1 of the SAUIS small is model stays inside the biologically feasible region Ω when β − δ − (β − βa )a∗0 − (β − βu )u∗0 < 0. (3) perturbed by ε, the three previous expressions must be positive. This occurs when β < δ and αa < δa . In any We point out that the expression on the left-hand side of other case, the SAUIS-ε model has no endemic equilibria the previous inequality is the same as the expression of bifurcating from P1 for ε > 0 small. In fact, the eigen- the unique eigenvalue of P2 for the unperturbed SAUIS values of P1 for the unperturbed SAUIS model are model that may take positive values (see equation (13) in [17] and comments surrounding). The other two eigen- λ1 (P1 ) = αa − δa , λ2 (P1 ) = −δu , λ3 (P1 ) = β − δ. values are either negative or have negative real part. In particular, this expression is negative if β < δ (since In consequence, P1 bifurcates to an endemic equilibrium ∗ β > βa , βu ), meaning that di is positive and, so, the for the system SAUIS-ε when P1 is hyperbolic stable in dε ε=0 the SAUIS model. Moreover, by the hyperbolic prop- SAUIS-ε system has an endemic equilibrium bifurcating erty the equilibrium remains stable for ε > 0 small in from P2 for ε > 0 small. When β > δ, the expression the SAUIS-ε model. The basic reproduction numbers can be positive or negative depending on the other pa- R0 = β/δ and R0a := αa /δa (see [17]) provide a clear rameters. This change can be controlled by taking βa as interpretation of this fact: for the SAUIS-ε model with a bifurcation parameter and so the bifurcation value is ε > 0 small enough, a stable equilibrium with all co- 1 ordinates positive and small emerges when R0 < 1 and βac := β − (β − δ − (β − βu )u∗0 ), a∗0 R0a < 1, corresponding to a non-spreading, dying out SAUIS epidemic. In this case, the small equilibrium val- as shown in [17]. In the SAUIS model, the system shows ues of aware, unwilling and infected are essentially fed by a transcritical bifurcation as βa passes through the bifur- 6 cation value and the authors illustrate that this bifurca- tion may occur in two different directions. That is, by changing the stability of P2 , from an stable equilibrium P2 for βa < βac a forward stable endemic equilibrium may bifurcate; or from an unstable equilibrium P2 for βa > βac a backward unstable endemic equilibrium may bifurcate. In both cases, the disease-free equilibrium is stable if βa < βac and unstable otherwise. In consequence, an endemic equilibrium bifurcates from the disease-free equilibrium in the SAUIS-ε when P2 is hyperbolic stable. As before in the case of P1 , the stability is preserved for ε > 0 small because of the hyperbolic property. The fact that the equilibria enter the region Ω for ε > 0 when they are hyperbolic stable is not surprising. Indeed, the vector field of the SAUIS-ε model at the boundary of Ω points towards its interior for ε > 0. This would be in contradiction with a hyperbolic unstable equilibrium entering Ω from the boundary. FIG. 1. Hopf-bifurcation curves in the (βa , αi ) parameter A similar treatment of the effect of immigration of in- space for ε = 0 (dashed line), ε = 10−5 (dotted line) and ε = fected individuals on the disease-free equilibrium of a 10−4 (solid line). Parameters: δ = 1, δa = 0.01, δu = 0.05, general epidemic model (without awareness), in terms β = 3, βu = 0.5, αa = 0.01, νa = 1. of the basic reproduction number, is given in [37]. Concerning the endemic equilibria P3 , by means of the implicit function theorem, we know that the root e∗ (ε) of (2) will persist inside Ω for ε > 0 small enough under the classical transversal condition as well as its stability. We refer to [38] for further information on the dynamical techniques used in this section and the forthcoming one. B. Robustness of the oscillatory regime The stability of an endemic equilibrium can change under a suitable election of parameters’ values. In par- ticular, for ε = 0, a Hopf-bifurcation curve H0 in the (βa , αi ) parameter space was obtained in [17]. Here we also do the analysis for ε = 10−5 and ε = 10−4 , which are small but still large enough to allow the existence of Hopf-bifurcation curves Hε clearly separated from H0 FIG. 2. Stability modulus of the Jacobian matrix at the en- (see Fig. 1). For each value of ε, the regime of sus- demic equilibrium of (2) as a function of αi . Parameters: tained oscillatory solutions of system (2) lies inside the δ = 1, δa = 0.01, δu = 0.05, β = 3, βa = 0.1, βu = 0.5, region of the parameter space limited by βa = 0 and αa = 0.01, νa = 1, ε = 10−4 . the corresponding Hopf-bifurcation curve. Outside this region, solutions tend to a stable endemic equilibrium. This behaviour is due to the fact that the real eigen- at the endemic equilibrium has a pair of pure imaginary value λr of the Jacobian matrix J of system (2) at the eigenvalues. Precisely, this condition is endemic equilibrium is always negative for any point of the considered region. Actually, λr is smaller than the c0 − c1 c2 = 0 with c1 > 0, (4) real part of the conjugate pair of complex eigenvalues λ± that constitute, together with λr , the spectrum σ of J. where c0 = − det(J), c1 is the sum of the principal minors So, this means that the stability modulus of J, namely, of J, and c2 = −trace(J). max{Re(λ) | λ ∈ σ(J)}, is given by Re(λ± ). From Fig. 1 we see that, for ε > 0, the region lim- We recall that to compute the Hopf-bifurcation curve ited by Hε is entirely contained in that limited by H0 . in the (βa , αi ) parameter space we need to find the solu- It follows that the more ε increases, the more the up- tions (βa∗ , αi∗ , a∗ , u∗ , i∗ ) of the system of equations given per part of the original region of the oscillatory regime by the three equilibrium equations of system (2) together is reduced, whereas the lower part of this region remains with the condition that follows from Theorem 2.1 and Ta- almost the same for the three bifurcation curves. This ble 1 in [39] which guarantees that the Jacobian matrix J fact suggests the existence of an abrupt transition from 7 the non-oscillatory regime below the lower branch of the parameters, p1 and p2 , running 50 experiments for each curves to the oscillatory one in the upper side of this pair (p1 , p2 ), when in addition the number of pairs is of branch. This transition is hardly perturbed by the exter- the order of 103 . In the particular case of regular ran- nal infection rates ε, provided that they are small enough. dom networks, this serious drawback can be overcome In fact, the qualitative behaviour of the solutions of by using what we will call the fast Gillespie algorithm system (2) presents significant differences between pa- (FGA in what follows). The FGA crucially relies on the rameter values near the upper and the lower boundary of following mean-field hypothesis (MFH): on a regular ran- the region determined by the Hopf-bifurcation curve. Al- dom network of big enough degree, the probability that though all solutions outside the confined zone are foci, a a neighbour of a node has a given state can be approx- simple study of the eigenvalues of the Jacobian matrix at imated by the fraction of nodes on the entire network the endemic equilibrium shows that the transition from having that particular state. Let us see how the MFH focus state (damped oscillations) to a periodic regime is can be used to speed up the Gillespie algorithm. faster through the lower boundary than through the up- Assume that during an experiment over a given net- per one. work a susceptible node n gets infected. In this case, More precisely, it follows that the absolute value of the Gillespie algorithm updates the state of n (from S the derivative of the stability modulus that controls the to I) and then explores all neighbours of n in order to oscillatory motion is much larger for lower values of αi update the number and type of links to be considered on the bifurcation curve (lower branch) than for higher in the next time step. For instance, if a neighbour of n ones (upper branch). In Fig. 2 we represent the stability is aware, we lose a link of type A − S with associated modulus for βa = 0.1 in terms of αi . This fact becomes weight αa0 + νa0 , and we gain a link of type I − S with crucial to understand the difference of sensibility between associated weight β 0 + αi0 . This exploration of the neigh- the two boundaries in the stochastic Hopf-bifurcation di- bours of a node through an adjacency matrix (usually a agrams appearing in next sections. A good example of pointer of pointers), that has to be done at every discrete this fact is the diagram presented in Fig. 9, where the time step, is one of the main computational loads of the simulations delimit the lower boundary of the Hopf curve classical Gillespie algorithm on graphs. almost identically as the theoretical curve. However, the But assume now that the approximation given by the upper boundary is hardly identified without taking into MFH assumption is good enough. Let N , NS , and NI account the amplitude of the signal. be respectively the total number of nodes, susceptible nodes, and infected nodes in a regular network of degree k. If the MFH holds, then the total number of links of IV. STOCHASTIC SIMULATIONS type I − S can be simply computed as kNI NS /(N − 1). Since the infection event I + S −→ I + I has rate β 0 = A. General simulation setup β/k, the total weight associated to all such events is then βNI NS /(N −1). Analogously, the total weight associated As usual in the setting of continuous-time stochastic to the event A+S −→ A+A, with rate αa0 = αa /k, would simulations, we use the well-known Gillespie algorithm on be αa NA NS /(N − 1), and so on. graphs [40]. All networks are randomly generated using The high speed of the FGA is achieved because the pro- the configuration model algorithm [41]. Given a network gram does not manage any particular network but only of size N , a combination of model parameters, and an ini- three integer variables NA , NI , NU (absolute numbers of tial condition a(0), u(0), i(0), we run 50 independent sim- aware, infected, and unwilling nodes respectively), with ulations, each corresponding to a random distribution of NS = N −NA −NI −NU . The total weight of all possible a(0)N , u(0)N , i(0)N , and (1−a(0)−u(0)−i(0))N nodes events is then having respectively the initial states of aware, unwilling, NI NS (β + αi ) NA NS (αa + νa ) NI NA βa infected, and susceptible. For any experiment, we store R := + + N −1 N −1 N −1 the evolution of a(t), u(t), and i(t) as three time series NI NU β u of 212 equally-spaced points in the interval [0, T ], where + + NI δ + NA δa + NU δu + (N − NI )ε. T is the maximum running continuous-time of the simu- N −1 lation. All along the paper, the caption of each reported A particular event as, for instance, I + S −→ I + I is figure obtained by simulation includes the specification chosen with probability (NI NS β/(N − 1))/R. In this of the values of N , a(0), u(0), i(0), and T . case, we just increase NI by 1, recompute R according to the previous formula and proceed to the next time step. There is no need to store a particular adjacency matrix B. The special case of regular random networks and explore the neighbours of any particular node, sim- ply because the MFH assumption allows us to work just It is worth noticing that the classical Gillespie algo- with the absolute numbers of aware, infected, and un- rithm on graphs is highly time-consuming when executed willing nodes. Observe that FGA is independent of the over a network of about 104 nodes, in such a way that it degree k. In other words, k is not a parameter of the is not feasible to construct a bifurcation diagram on two algorithm. The FGA performs statistically exact simula- 8 tions of Markovian epidemic processes over regular ran- We can observe this fact in Fig. 5, which shows that dom networks of high enough degree (that is, as long as all the region of the (βa , αi )-space where periodic solu- the MFH applies). It is worth mentioning that what tions (interior of the Hopf-bifurcation curve) and weakly we have called FGA can be identified with the original damped oscillatory solutions are predicted by system (2) version of the algorithm [40], which was aimed at the lies within the extinction zone (dark region). Note that, stochastic simulation of a fully mixed chemically react- for a fixed value of αi , the greater βa is, the higher the ing system. disease prevalence at the endemic equilibrium and, hence, Let us see to which extent the outputs of the FGA and the lower the extinction probability of the epidemic. Con- the Gillespie algorithm are essentially equivalent. Recall versely, for a fixed βa , the higher the alerting rate αi is, (Subsection III B) that for δ = 1, δa = 0.01, δu = 0.05, the lower the prevalence because aware individuals are β = 3, βu = 0.5, αa = 0.01, νa = 1, and ε = 10−4 , more easily created and, hence, the higher the extinction a Hopf-bifurcation curve was obtained in the (βa , αi )- probability is. This is the reason why the non-extinction plane. In Fig. 3 we have shown the time evolution of region corresponds to the lower right part of this figure. the fraction of infected nodes according to the numeri- Fig. 6a shows one stochastic simulation of an oscil- cal integration of system (2), together with the adaptive lating epidemic on a random regular network of 1000 averaged outputs (see Section IV D) of the Gillespie al- nodes and, also, the moments where external infections gorithm over a regular random network of size N = 1000 are introduced. Such an oscillatory behaviour is, in fact, and degree 50, and FGA with N = 1000. When αi = 0.04 predicted by system (2). However, the extremely low and βa = 0.1 (Fig. 3a), we are below the Hopf curve and prevalence attained in each cycle makes the stochastic so we have an stable endemic equilibrium. Observe that extinction of the disease unavoidable, unless the occa- the three curves are essentially identical. For αi = 0.14 sional introduction of imported cases takes place when and βa = 0.1 (Fig. 3b) we are inside the Hopf curve and population awareness is low. The crucial role of such we get a stable periodic orbit. In this case, the outputs an introduction in maintaining the epidemic is revealed from Gillespie and FGA seem qualitatively equivalent up in Fig. 6b, where the arrival of imported cases has been to stochastic fluctuations. forced to cease (ε = 0) just at the beginning of the fourth Of course, we should give a precise meaning to the sen- flare-up. As expected, the fifth flare-up (dashed line) that tence seem qualitatively equivalent up to stochastic fluc- would appear by keeping ε = 10−4 now does not happen tuations. This is in fact the aim of Section IV D, where and the epidemic dies out. The total number of imported we give a detailed explanation about the statistical treat- cases from t = 0 to t = 1200 is 119 in Fig. 6a, whereas it ment of the simulation data in order to test the signif- is equal to 85 in Fig. 6b where ε = 0. icance of the oscillatory regime. In Fig. 4 we show the The low number of imported cases in the previous ex- complete Hopf-bifurcation diagram for the detection of ample, about 1 case every 10 infectious periods on av- the oscillatory regime in (βa , αi ) ∈ [0, 0.7] × [0, 1] after erage in a population of size 1000, shows that it is not processing the output data obtained by both the classi- necessary to have a high number of imported cases to cal Gillespie algorithm over a regular random network of prompt the occurrence of important flare-ups in popu- N = 1000 nodes and degree 50 (Fig. 4a) and the FGA lations whose individuals have a low level of awareness. with N = 1000 (Fig. 4b). We stress that this figure is This situation reminds, for instance, of what happened intended only to compare both algorithms. In particular, in New Zealand after the praised management of their we have chosen here N = 1000 since a higher order for N first wave of COVID-19, where a reemergence of cases makes the computation of the Hopf-bifurcation diagrams occurred after weeks with no community cases once so- under the Gillespie algorithm on networks highly costly cial restrictions were lifted [42]. in time. Observe that the two diagrams are essentially identical, showing that, as expected, the FGA is a good substitute of the Gillespie algorithm on regular graphs D. Statistical significance of the oscillatory regime even for degrees as small as 50 (over 1000 nodes). Our approach for the statistical detection of oscilla- tions in the time series has two parts. The first one con- C. Impact of imported cases in the prevention of sists in the generation of three signals from the output stochastic extinctions of the stochastic simulations, one for each state variable of the model. As stated in Section IV A, for each com- Imported cases have an impact on the evolution of an bination of parameters we run K = 50 independent ex- epidemic if they happen when both disease prevalence periments. Each realisation consists on 212 points, of and awareness level in the population are low. In their which we eliminate the first half seeking for stationarity absence (ε = 0), an epidemic becomes eventually extinct of the time series, ending up with M = 211 points. In because the number of infected nodes becomes extremely stochastic epidemic models, the usual way to construct a low when the level of awareness in the population is very single signal from data is through the average. It is well high. In such a situation, it is well known that stochastic known that ODE systems are a good approximation of extinctions are extremely likely. the mean of realisations of stochastic processes in systems 9 (a) (b) FIG. 3. Evolution of i(t) for t ∈ [1500, 3000] predicted by system (2) and outputs from Gillespie and FGA (adaptive mean over 50 experiments). Parameters of the simulations: N = 1000, degree k = 50 for the Gillespie algorithm, δ = 1, δa = 0.01, δu = 0.05, β = 3, βu = 0.5, αa = 0.01, νa = 1, ε = 10−4 , βa = 0.1, i(0) = 0.1, a(0) = u(0) = 0.2. (a) αi = 0.04, (b) αi = 0.14. (a) (b) FIG. 4. Hopf diagram of SAUIS-ε obtained according to the description in Section IV D, using (a) Gillespie algorithm and (b) FGA with N = 1000 nodes and time T = 3000. Parameters: degree k = 50 (for the Gillespie algorithm), δ = 1, δa = 0.01, δu = 0.05, β = 3, βu = 0.5, αa = 0.01, νa = 1, ε = 10−4 , i(0) = 0.1, a(0) = u(0) = 0.2, 50 experiments for each pair (βa , αi ). The gradient of colours evidences the amplitude of the signal corresponding to the fraction of infected nodes. The black line is the theoretical Hopf-bifurcation curve. with a large number of components (nodes, molecules in measures against infection. Such periods are followed by reaction systems, etc.). However, although the standard relaxation in the adoption of preventive measures because average of trajectories works for endemic equilibria with of a decay in awareness which, in turn, is translated into a high prevalence (see Fig. 3a), it is not always a good new epidemic flare-ups. During these time intervals of option when a system exhibits dynamics where some of low prevalence (see Fig. 6a), the stochastic trajectories of its components are present in low numbers [43]. the epidemic show significant random fluctuations which In particular, the relationship between behaviour and lead to a phase shift in the time series with respect to disease we are considering allows the existence of epi- the deterministic trajectory given by the solution of the demics with periods of very low prevalence once individ- ODE model. Consequently, when the average of trajec- uals have become aware and adopt efficient protective tories is computed, the resulting signal has an important 10 the resulting data Y . To do so we define the periodogram of the zero-mean time series Ŷ := Y − Y by 2 M ∆t X P (ω0 ) = Ŷ (j) , M j=1 2 M 2∆t X 2πni P (ωn ) = Ŷ (j)e− M (j−1) M j=1 for n = 1, . . . , M 2 − 1, and 2 M ∆t X P (ωM/2 ) = Ŷ (j)e−πi(j−1) , M j=1 FIG. 5. Fraction p of stochastic simulations with positive with ωn = Mn∆t and ∆t being the time-step of the time prevalence up to time T = 3000 for ε = 0 and N = 10000 us- series. The periodogram is suitable to find hidden period- ing FGA algorithm. The Hop-bifurcation curve (white curve icities on the data. Indeed, data following a pure random within the extinction region) is included for a better visuali- process show a flat periodogram. However, if the time sation of the extinction range. Parameters: δ = 1, δa = 0.01, series presents a periodic motion then the frequencies in- δu = 0.05, β = 3, βu = 0.5, αa = 0.01, νa = 1, ε = 0, volved are shown as peaks of the periodogram. Fisher’s i(0) = 0.1, a(0) = u(0) = 0.2, 100 experiments for each pair (βa , αi ). g-test [44] is commonly used in this scenario to reject the null hypothesis that the random process is Gaussian white noise against the alternative hypothesis that the decrease of the amplitude of the oscillations and even a series contains a deterministic periodic component of un- deformation of the periodic component. specified frequency. The statistic taken into account to To overcome this difficulty, we perform an adaptive reject the null hypothesis is the maximum of the peri- mean of the signals. The idea is to locally align the sig- odogram compared with the sum of all the values of the nals before computing the average, which acts in favour of periodogram. That is, preserving periodic motion when the mean is performed. maxn P (ωn ) The first part of this method is classic in signal process- g = PM/2 . ing. To begin with, we perform a moving average com- n=1 P (ωn ) puted over a sliding window of length 21 centred about each element in the time series. This acts as a low-pass The previous maximum, as well as the summation, is filter to each time series attenuating the signal and omit- taken for n = 1, . . . , M/2. This is so because the term ting extreme frequencies. Then, taking the first time P (ω0 ) contains only information about the mean of the series, say X1 , as a sample, we determine the time delay data. This term plays no role in the detection of oscil- between X1 and each of the other signals. To accomplish latory motion and, in fact, vanishes for zero-mean time this, for each time series Xk , k = 1, . . . , K, we compute series. The test is performed by computing the realised the cross-correlation between X1 and Xk . The position value g ∗ of g from the data and then use the exact dis- τk where the maximum is reached, that is tribution of g given by Fisher in [44], b  M X (−1)k−1 M 2 −1 τk := argmaxn ((X1 ? Xk )(n)), n = 1, . . . , M, 2 !(1 − kx) P (g > x) = M , k=1 k! 2 − k ! corresponds to the position where the signals are best aligned. For each of the signals Xk we define a new one, where b is the largest integer less than 1/x. If the prob- X̂k , by applying a circular shift of τk positions to Xk . In ability P (g > g ∗ ) is less than α then the null hypothesis practice the previous cross-correlation is computed for can be rejected at level α. For further details on Fisher’s small delays to ensure that the alignment is local. The test for hidden frequencies we refer to [45–47]. adaptive mean we consider is given by For the time series resulting from the parameters we K are interested in, the direct application of the Fisher’s 1 X test results positive for all configurations. This is so for Y := X̂k . K two reasons. First, all motions of the ODE system (2) k=1 have an oscillatory component since they correspond to The second part of the detection of oscillations consists either a focus or a periodic orbit. Second, for the expres- in proving the statistical significance of the periodicity of sion of the probability above, it easily follows that the 11 (a) (b) FIG. 6. (a) Time evolution of a stochastic simulation of an epidemic on a random regular network of size N = 1000 and degree k = 100 showing the input of imported cases (infections from abroad) for ε = 10−4 (black dots on the time axis). (b) Fraction of infected nodes in the same simulation until the beginning of the fourth flare-up where ε = 0 (lockdown). Dashed line: fraction of infected nodes without lockdown. Parameters: δ = 1, δa = 0.01, δu = 0.05, β = 3, βa = 0.05, βu = 0.5, αi = 0.15, αa = 0.01, νa = 1, i(0) = 0.1, a(0) = u(0) = 0.2. larger the number of points M , the smaller the value of namely, 0.01 (and 0.03 when N = 1000 as in Figs. 4 the probability. To be more precise in the determination and 10). This value comes from the so-called linear noise of oscillatory motions, we reduce the points of the signals approximation used in modelling of chemical reaction ki- after the adaptive mean in a ratio 1 : 24 . That is, we end netics [48, 49]. Such an approximation assumes that, in up with M = 27 points. This procedure enables to detect a system formed by different “chemical species” (nodes the strongest periodic motions in the parameter space, in different states in our context), the standard devia- since weaker periodic signals will not pass Fisher’s test tion of random fluctuations about the mean number of with fewer points. Fig. 7 illustrates it with two param- molecules of these chemical species scales as the square eter configurations. On the left, periodograms with 211 root of the size of the system (the number N of nodes in points pass Fisher’s test for both configurations. On the our network). So, dividing these mean numbers by N , right, the periodograms of the same averaged signals with it follows that, under this approximation, the standard 27 points. In this case the configuration corresponding to deviation of the random fluctuations of the fractions of the periodogram on the top side does not pass Fisher’s A, U , and I nodes is proportional to N −1/2 . test, whereas the one on the bottom side still passes the test. Notice that these parameter configurations corre- spond to the signals of top panels in Fig. 11. It is worth E. Oscillations in other network architectures to mention that the frequencies avoided by the reduc- tion of points are much larger than the frequencies of the Another aspect of the robustness of the oscillations is data as can be seen in Fig. 7 and information about the their likelihood when other network topologies are con- periodicity of the signal is not lost. sidered. Are they still present when the epidemic spreads For those parameters passing the Fisher’s g-test with on more heterogeneous networks? If this is the case, are a p-value less than or equal to α = 0.01 we consider they observed in the solutions of system (1)? To answer the corresponding peak frequency of the time series. The these questions, we have used the configuration model al- bifurcation diagram shows those parameters with an esti- gorithm to generate a Poisson network with mean degree mated frequency larger than the minimum observable fre- 50 and an exponential network with mean degree 50 and quency fmin = 1/1500. A gradient of colours illustrates minimum degree 25. The first degree distribution corre- the amplitude of the averaged signal corresponding to sponds to the well-known Erdös-Rényi random graphs, infected nodes. In Fig. 9 we have shown the correspond- whereas exponential degree distributions have been ob- ing bifurcation diagram for the outputs of the SAUIS-ε served in empirical contact networks [50] and have a model using the FGA with N = 10000. The thresh- much higher variability. old below which oscillations are considered to be inter- Unfortunately, in this case it is not possible to reduce nal random fluctuations (noise) of the system is given by (1) to a simpler system, as it was done in Subsection III B the inverse of the square root of the number of nodes, for regular random networks. So, the algebraic approach 12 (a) (b) (c) (d) FIG. 7. Periodograms of SAUIS-ε obtained using FGA with N = 10000 nodes and time T = 3000 (adaptive mean over 50 experiments). The number of points used is M = 211 for (a) and (c), and M = 27 for (b) and (d), according to the procedure in Section IV D. Parameters: δ = 1, δa = 0.01, δu = 0.05, β = 3, βa = 0.1, βu = 0.5, αa = 0.01, νa = 1, ε = 10−4 , i(0) = 0.1, a(0) = u(0) = 0.2. (a) and (b): αi = 0.04. (c) and (d): αi = 0.14. used there for computing Hopf-bifurcation curves is no pute CVx(t) = σx(t) /x̄ × 100 with longer feasible. Instead, using βa and αi as tuning pa- Z T2 rameters (with increments of size 0.002 or even smaller 1 when approaching the turning point of the curve), we x̄ = x(t) dx T2 − T1 T1 have obtained an approximation to these curves by nu- merically integrating system (1) for N = 1000 (i.e. the and full system of 3000 equations) until a long enough time Z T2 T2 . Then, for each pair (βa , αi ), we computed the coeffi- 2 1 2 σx(t) = (x(t) − x̄) dx. cient of variation (CV ) of the fraction of aware, unwilling T2 − T1 T1 and infected nodes from the solution Ai (t), U i (t), and I i (t) of (1) for the last 1000 units of time (T P2 −i T1 = Finally, an approximate Hopf-bifurcation curve is ob- 1000). Precisely, for the fraction x(t) = i x (t)/N tained by using the criterion that, for a given pair where xi (t) = Ai (t), U i (t), I i (t), respectively, we com- (βa , αi ), a solution is not periodic if the corresponding CV ≤ 0.1% for the three fractions of nodes. Note that, close to the Hopf-bifurcation, solutions classified as peri- odic when we integrate the system up to a given time can 13 gebraic Hopf-bifurcation curve, do not depend on the de- gree, while system (1) does depend on it. Fig. 8 shows that, as expected, increasing the integration time from T2 = 10000 to T2 = 20000 leads to a better approxi- mation to the whole algebraic curve (solid line). This figure also shows that, for T2 = 20000, the disagreement with respect to the algebraic curve is only perceptible in its upper branch. Therefore, the (approximate) Hopf- bifurcation curves for Poisson and exponential networks were constructed from the solutions of system (1) using an integration time T2 = 20000 and a threshold value of CV equal to 0.1% for the fraction of the three types of non-susceptible nodes. F. Simulation results and discussion FIG. 8. Comparison of the algebraic Hopf-bifurcation curve for regular random networks with two approximations ob- The values of the parameters used in the simulations tained integrating system (1) up to two different times and are taken from [17] (except for the rate of imported cases using for CV = 0.1% to discriminate between periodic and ε) and reflect what we consider it is a natural scenario, non-periodic solutions. although they are not intended to model any particu- lar disease. First, aware (and unwilling) individuals are affected by lower transmission rates because of the adop- tion of preventive measures. Second, the mean infec- tious period is much shorter than the mean duration of awareness. Finally, it is assumed that is more difficult for an aware individual to convince a susceptible one to become aware than to convince him to become simply unwilling, i.e., to adopt preventive measures but with- out willingness to convince others about the risk of in- fection. The values of the awareness and unwillingness decay rates lead to oscillations whose period is about 1.3 times the sum of the mean duration of the awareness period (1/δa = 100) and the mean duration of the un- willingness period (1/δu = 20). The procedure of the adaptive mean of the signals from the stochastic simulations gives an excellent reconstruc- tion of the oscillatory behaviour of single trajectories. In Fig. 11, we show the comparison of a solution of system FIG. 9. Hopf diagram obtained according to the description (2), the adaptive mean of the 50 realisations of an epi- in Section IV D of the FGA with N = 10000 nodes and time demic on a regular random network of size N = 10000 T = 3000. The gradient of colours evidences the amplitude generated by means of the FGA, and the trajectory used of the averaged signal of the fraction of infected nodes. The as the reference for the alignment of the signals. Note black line is the algebraic Hopf-bifurcation curve for ε = 10−4 . that this alignment preserves the periodicity of the origi- nal signals (trajectories) but, at the same time, it reduces their amplitudes. With this respect, it is worth recalling become non-periodic when longer times are considered. that the amplitudes shown in the Hopf diagrams do not This is particularly relevant when trying to delimit the correspond to stochastic simulations of single epidemics, upper branch of the bifurcation curve. Here the transi- but to the adaptive means of 50 stochastic trajectories tion from periodic solutions to very weakly damped solu- for each pair (βa , αi ) of parameter values (see Sect. IV D tions is hardly noticed because of its flatness, which was for details). already observed when computing the algebraic curve for Hopf-bifurcation diagrams show that the abrupt tran- the regular case (see Fig. 2). sition between a non-oscillatory regime and an oscillatory To calibrate our criterion with respect to the integra- one given by the lower branch of the Hopf-bifurcation tion time, we compared the algebraic Hopf-bifurcation curve is quite accurately captured in all the settings we curve for system (2) with those obtained from system have considered. The robustness of this transition is, in (1) for a regular random network of degree 50. No- fact, already observed when this curve is computed for tice that system (2) and, hence, the corresponding al- different rates of imported cases (see Fig. 1). 14 (a) (b) FIG. 10. Hopf diagram of SAUIS-ε obtained according to the description in Section IV D with N = 1000 nodes and time T = 3000, in two different network architectures of mean degree 50: (a) Poisson and (b) Exponential. Parameters: δ = 1, δa = 0.01, δu = 0.05, β = 3, βu = 0.5, αa = 0.01, νa = 1, ε = 10−4 , i(0) = 0.1, a(0) = u(0) = 0.2, 50 experiments for each pair (βa , αi ). The gradient of colours evidences the amplitude of the signal corresponding to the fraction of infected nodes. The black line is the theoretical Hopf-bifurcation curve. The best agreement is achieved for regular random which has been claimed to describe contact patterns in networks of size N = 10000. Fig. 9 shows that the re- real populations. Interestingly, in all cases the predicted gion with periodic solutions with amplitudes larger than region of the parameter space for the oscillatory regime 0.05 falls neatly within the limits predicted by the Hopf- satisfy that βa < βu , a natural condition according to bifurcation curve. Comparing this figure with the panels the higher level of alertness assumed for aware individu- in Fig. 4 obtained for N = 1000, one can see that the als with respect to the unwilling ones. Moreover, these agreement of the observed region of periodic solutions regions always contain those averaged stochastic trajec- with high amplitudes (yellow and orange squares) with tories with the highest amplitudes. the predicted one clearly increases with the number of nodes. Moreover, we recall that all the simulations have been carried out until time T = 3000 which can be not V. CONCLUSIONS enough for weakly damped oscillations to disappear. A similar agreement with the predictions as the one in The interplay between human behaviour and epidemic Fig. 4 is also observed in networks with different topolo- spread has been considered in many papers during the gies but with the same number of nodes, N = 1000 (see last 15 years [2, 51]. Most of them have as the main Fig. 10). Moreover, we can observe that the differences goal to elucidate the dependence of the basic reproduc- between the Hopf-bifurcation curves when the degree dis- tion number on the different behavioural responses as, tributions have different variability are remarkable. For for instance, social distancing, rewiring of connections, regular random networks and Poisson networks, these dif- awareness, etc. [10, 11, 13, 14, 23, 25]. Few of them ad- ferences are small, i.e. both types of networks lead to dress the existence of oscillating epidemics [20, 21, 27–29] similar (but not equal) Hopf-bifurcation curves. How- and, as far as we know, only in [17] sustained oscillations ever, the resulting curve for an exponential network is arise uniquely from the interaction between awareness clearly different, with a turning point at a higher value dissemination and epidemic spread, i.e. without the need of αi and lower values of βa (see Fig. 10b). Interestingly, of any recruitment of new (susceptible) individuals. the region of the (βa , αi ) parameter space correspond- In this paper we have challenged the existence of oscil- ing to periodic solutions with large amplitudes (yellow lations predicted by the SAUIS model introduced in [17] squares and most of the orange ones) clearly falls inside by means of stochastic simulations on different network the region limited by the Hopf-bifurcation curve. topologies. The model has been formulated on networks All in all, stochastic simulations confirm the predic- and considers the existence of a very small inflow of im- tions of the ODE model showing that the proposed mech- ported cases. These cases are essential to keep the epi- anism for prompting awareness is able to generate recur- demic going on because disease prevalence attains very rent epidemic cycles on different network architectures, low levels when there is a high degree of awareness in an including the one with an exponential degree distribution oscillating epidemic. This is not a problem at all in a 15 (a) (b) (c) (d) FIG. 11. Evolution of i(t) for t ∈ [1500, 3000] predicted by system (2) (solid line), adaptive mean of outputs from FGA (over 50 experiments, dashed line) and a single experiment from FGA (dotted line). Parameters of the simulations: N = 10000, δ = 1, δa = 0.01, δu = 0.05, β = 3, βa = 0.1, βu = 0.5, αa = 0.01, νa = 1, ε = 10−4 , i(0) = 0.1, a(0) = u(0) = 0.2. (a) αi = 0.04, (b) αi = 0.14, (c) αi = 0.16, (d) αi = 0.26. deterministic framework, but it leads to an unavoidable creases (N = 10000), periodic epidemics with amplitudes stochastic epidemic extinction in all parameter combina- larger than 5% carried out on regular random networks tions where sustained oscillations are present. fall within the oscillatory regime in the parameter space predicted by the ODE model (see Fig. 9). So, the robust- Our simulations show that the presence of a small num- ness of the predictions about the existence of oscillations ber of imported cases allows, when the disease prevalence due to pure behavioural changes has been established. is low and the number of aware and unwilling individu- als decreases, an oscillatory behaviour of the stochastic An interesting example of oscillating epidemics is given epidemics which closely resembles the one predicted by by the evolution of the incidence rate of sexual transmit- the deterministic counterpart of the model. In particular, ted diseases (STDs) during the last decades. The reemer- the existence of an abrupt transition from a stationary gence of STDs like syphilis and gonorrhea occurring since regime where oscillations are strongly damped to an os- the mid-1990s [52] has been associated with a decrease cillatory one with amplitudes of more than 10% in the in awareness after the introduction of the antiretroviral number of infected nodes predicted by the ODE model therapy for HIV and, indeed, it appeared after an inci- is clearly observed in the simulations over different types dence decline in the 1980s. This decline coincided with of networks. Moreover, when the size of the network in- the emergence of the global AIDS pandemic and has been 16 attributed to preventive behavioural changes in response ACKNOWLEDGMENTS to HIV campaigns during that time [53]. However, it has been also claimed that the drop in incidence before 1984 D.J. and D.R. have been partially supported by the occurred too early to be ascribed to such induced be- Agencia Estatal de Investigación and Ministerio de Cien- havioural changes and may be part of a long term cyclic cia, Innovación y Universidades grant MTM2017-86795- trend of this type of diseases (although the nature of this C3-1-P, and D.J. and J.S. have been partially supported periodic behaviour remains unexplained) [54]. by the grant PID2019-104437GB-I00 of the Ministerio de Ciencia e Innovación of the Spanish government. D.J. and D.R. are members of the Consolidated Research Group 2017 SGR 1617 funded by the Generalitat de Catalunya. J.S. is member of the Consolidated Research Group 2017 SGR 01392 of the Generalitat de Catalunya. [1] N. Ferguson, Nature 446, 733 (2007). [23] C. Llensa, D. Juher, and J. Saldaña, J. Theor. Biol. 352, [2] S. Funk, M. Salath, and V. Jansen, J. R. Soc. Interface 71 (2014). 7, 1247 (2010). [24] S. Risau-Gusmán and D. Zanette, J. Theor. Biol. 257, [3] Centers for disease control and prevention. Sexually 52 (2009). transmitted disease surveillance 2018. Atlanta: U.S. De- [25] D. Zanette and S. Risau-Gusmán, J. Biol. Phys. 34, 135 partment of Health and Human Services (2019), Avail- (2008). able at https://www.cdc.gov/std/stats/. [26] S. V. Scarpino, A. Allard, and L. Hébert-Dufresne, Na- [4] Global measles and rubella strategic plan: 2012–2020. ture Physics 12, 1042 (2016). world health organization. ISBN: 978 92 4 150339 6. [27] A. Szabó-Solticzky, L. Berthouze, I. Kiss, and P. Simon, (2012). J. Math. Biol. 72, 1153 (2016). [5] B. Renner, M. Gamp, R. Schmälzle, and H. Schupp, [28] C. Bauch, Proc. R. Soc. 272, 1669 (2005). Health risk perception, in International Encyclopedia of [29] J. Velasco-Hernández, F. Brauer, and C. Castillo-Chaves, the Social & Behavioral Sciences (Oxford: Elsevier, New IMA J. Math. Appl. Med. Biol. 13, 175 (1996). York, 2015) pp. 702–709, 2nd ed. [30] Measles in England. Blog: Public health [6] V. Capasso and G. Serio, Math. Biosc. 42, 43 (1978). matters, Public Health England (2019), [7] W. Liu, H. Hethcote, and S. Levin, J. Math. Biol. 25, https://publichealthmatters.blog.gov.uk/2019/08/19/ 359 (1987). measles-in-england/. [8] X. Lu, S. Wang, S. Liu, and J. Li, Math. Biosci. Eng. 14, [31] Measles: Four european nations lose eradi- 1317 (2017). cation status. BBC News (29 August 2019), [9] S. Ruan and W. Wang, J. Differential Equations 188, https://www.bbc.com/news/health-49507253. 135 (2003). [32] Url: https://www.statista.com/statistics/1106879/number- [10] S. Funk, E. Gilad, and V. Jansen, J. Theor. Biol. 264, of-coronavirus-cases-in-iceland-by-origin-of-infection/ 501 (2010). (2020). [11] I. Kiss, J. Cassell, M. Recker, and P. Simon, Math. Biosci. [33] Url: https://www.smh.com.au/world/asia/singapore- 225, 1 (2010). sees-huge-surge-in-new-virus-cases-20200419- [12] C. Granell, S. Gómez, and A. Arenas, Phys. Rev. Lett. p54l53.html (2020). 111, 128701 (2013). [34] S. Gómez, A. Arenas, J. Borge-Holthoefer, S. Meloni, and [13] D. Juher, I. Kiss, and J. Saldaña, J. Theor. Biol. 365, Y. Moreno, Europhys. Lett. 89, 38009 (2010). 457 (2015). [35] P. V. Mieghem, Computing 93, 147 (2011). [14] F. Sahneh, F. Chowdhury, and C. Scoglio, Scientific Re- [36] F. Sahneh, C. Scoglio, and P. V. Mieghem, IEEE/ACM ports 2, 632 (2012). Transactions on Networking 21, 1609 (2013). [15] F. Sahneh and C. Scoglio, in 50th IEEE Conference on [37] R. M. Almarashi and C. C. McCluskey, J. Math. Biol. Decision and Control (2011). 79, 1015 (2019). [16] S. Funk, E. Gilad, C. Watkins, and V. Jansen, PNAS 21, [38] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, 6872 (2009). Dynamical Systems, and Bifurcations of Vector Fields [17] W. Just, J. Saldaña, and Y. Xin, J. Math. Biol. 76, 1027 (Springer-Verlag, New York, 1983). (2018). [39] J. Guckenheimer, M. Myers, and B. Sturmfels, SIAM J. [18] B. M. Althouse and L. Hébert-Dufresne, J. R. Soc. In- Numer. Anal. 34, 1 (1997). terface 11, 20140575 (2014). [40] D. Gillespie, Annual Review of Physical Chemistry 58, [19] T. Britton, D. Juher, and J. Saldaña, Bull. Math. Biol. 35 (2007). 78, 2427 (2016). [41] T. Britton, M. Deijfen, and A. Martin-Löf, Journal of [20] T. Gross, C. D’Lima, and B. Blasius, Phys. Rev. Lett. Statistical Physics 124, 1377 (2006). 96, 208701 (2006). [42] M. Daalder, The NZ strains: Our second wave (5 Octo- [21] T. Gross and I. G. Kevrekidis, EPL 82, 38004 (2008). ber 2020), https://www.newsroom.co.nz/the-nz-strains- [22] D. Juher, J. Ripoll, and J. Saldaña, J. Math. Biol. 67, our-second-wave. 411 (2013). [43] S. Hahl and A. Kremling, Front. Genet. 7, 157 (2016). 17 [44] R. Fisher, Proc. R. Soc. A. 125, 54 (1929). [51] F. Verelst, L. Willem, and P. Beutels, J. R. Soc. Interface [45] P. Brockwell and R. Davis, Times Series: Theory and 13, 20160820 (2016). Methods, 2nd ed. (Springer-Verlag, New York, 1991). [52] D. Williamson and M. Y. Chen, N. Engl. J. Med. 382, [46] M. Priestley, Spectral analysis and time series (Academic 2023 (2020). Press, London, 1981). [53] K. Fenton and C. Lowndes, Sex. Transm. Infect. 80, 255 [47] A. Schuster, Terr. Mag. Atmos. Elect. 3, 13 (1898). (2004). [48] R. Grima, J. Chem. Phys. 133, 035101 (2010). [54] M. S. Green, E. Anis, D. Gandacu, and I. Grotto, Sex. [49] N. V. Kampen, Stochastic processes in Physics and Transm. Infect. 79, 116 (2003). Chemistry (Elsevier, Amsterdam, 2007). [50] S. Bansal, B.T.Grenfell, and L.A.Meyers, J. R. Soc. In- terface 4, 879 (2007).