Robustness of behaviorally induced oscillations in epidemic models under a low rate of imported cases

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 extinctions 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 simulations confirm the existence of sustained oscillations for different types of random networks, with a sharp transition from a nonoscillatory asymptotic regime to a periodic one as the alerting rate of susceptible individuals increases from very small values. This abrupt transition to periodic epidemics 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.


I. INTRODUCTION
The importance of the interplay between epidemic spreading and preventive behavioral responses in a globalized world has long been recognized and was specially highlighted after the SARS outbreak of 2003 [1,2]. The rise of the incidence rate of sexually transmitted diseases (STDs) [3] and the current resurgence of measles [4] are also examples of such an interplay. For STDs, increasing high risk sexual behavior and novel sexual networks are among factors responsible for their reemergence, whereas vaccine hesitance and distrust in public health intervention programs are among behavioral factors responsible for the rise of diseases like measles.
Risk perception is an important determinant of selfinitiated, voluntary protective behavior [5]. It constitutes the basic ingredient in many epidemic models to encapsulate human behavior in their formulation [2]. For instance, some extensions of classic deterministic compartmental models include the impact of behavior on disease transmission by assuming more general incidence rates than the standard one (bilinear). The latter is proportional to the product of the number of susceptible (S) and infectious individuals (I), βSI, whereas its generalizations assume a saturation with respect to the number of infectives to model a reduction of the contact rate in the presence of a high disease prevalence [6][7][8][9].
Other model extensions take into account awareness transmission among individuals. For example, one of them [10,11] divides each epidemic compartment [S, I, and R (recovered/immune)] into two subcompartments of aware and unaware individuals, respectively, and introduces the corresponding transition rates between subcompartments. * david.juher@udg.edu † david.rojas@udg.edu ‡ joan.saldana@udg.edu Other models assume one or more additional compartments consisting of aware (A) individuals [12][13][14][15]. Some of these works consider several awareness levels resulting from the assumption of a degradation in the quality of the information as it is passed from one individual to another [16,17]. In all these examples, the effect of preventive behavior is to modify the values of the epidemic parameters like the probability of infection or the recovery rate.
When the contact network structure of a population is explicitly considered, the effect of behavioral responses can also affect the contact structure itself when modeling social avoidance behaviors. This reduction of the exposure to disease has been modelled by means of preventive disconnection from infectious neighbors [18][19][20][21][22][23][24][25] or, also, by replacing some infected nodes by healthy ones [26], leading, in both cases, to dynamical networks. Here the assessment of disease prevalence is based on the individual neighborhood (local contacts), in contrast to homogeneous compartmental models where information about the prevalence is assumed to be globally available [2].
Under the previous modeling approaches, protective behavioral responses are triggered by the disease prevalence. As long as these responses are based on the global prevalence, one expects the likelihood of epidemic oscillations to be high. Moreover, the linear dependence of the standard incidence rates on the number of infected individuals implies that, when these oscillatory solutions occur, they should pass through low prevalence levels due to the lack of an abrupt switching behavior. A low prevalence, in turn, will drive the number of aware individuals down as a consequence of a lower perception of the contagion risk, and the cycle repeats again with a new rise in the number of infectives. In fact, several ODE epidemic models with transmission of awareness [17,27] or assuming self-initiated, voluntary vaccination [28] exhibit such periodic solutions under some values of the parameters. In particular, the variability in the propensity of aware individuals to further propagate awareness as the driving mechanism for sustained oscillations was proved in Ref. [17]. More precisely, it was proved there that, in the absence of demographics, standard incidence rates of infections and awareness do not lead to oscillations even though aware individuals disseminate awareness among susceptible ones, unless a class of individuals with a lower level of awareness (labeled U for unwilling to disseminate) is also assumed.
In this paper we study the robustness of the deterministic oscillations of a model (which we call SAUIS-ε) that is an extension of the SAUIS model studied in Ref. [17]. The reason for this choice is twofold. First, because the only driving mechanism for the existence of periodic solutions in the SAUIS model is the variability in the propensity of alerted individuals to further propagate awareness in the population. Second, because the model does not consider recruitment of susceptible individuals in terms of newborns or in terms of a prevalence-dependent recruitment of them into a core group, which turns out to be an essential requirement for having sustained oscillations in many epidemic models with and without vaccination [9,28,29].
Another behavioral mechanism also responsible for the occurrence of periodic solutions in epidemic models is preventive rewiring in adaptive contact networks. Such a mechanism aims to minimise the infection risk of susceptible individuals while maintaining them connected to the network [18,20,21,27]. So, both approaches try to explain epidemic oscillations on a pure behavioral basis.
In addition, damped epidemic oscillations have been predicted in Ref. [16], where an SIR epidemic model without demography is coupled with a sophisticated mechanism of degradation of the information quality. Such degradation is translated into individuals with different levels of awareness in the population. However, the eventual depletion of susceptible nodes prevents the occurrence of sustained oscillations.
In countries where an infectious disease has been declared 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 enough. These new cases are usually imported either from tourists, foreign workers, etc., or by local individuals that have been infected abroad in regions where the disease is endemic or large outbreaks are taking place. However, if the vaccination coverage in these countries decreases because of lower levels of awareness (for instance, in countries where vaccination is not mandatory), then such few cases can mount into major outbreaks.
In fact, regions that have achieved the WHO eradication status for a given disease can lose it as a consequence of a marked increase in the number of confirmed cases. For measles, a country attains this status when there is no endemic transmission for 12 months in a specific geographic area. For instance, the UK achieved WHO measles eradication status in 2017, based on data from 2014 to 2016. However, two years later, this status was lost after 991 confirmed cases in England and Wales in 2018, more than three times the number of cases (284) in 2017 [30]. Other European countries that have also lost the WHO eradication status for measles in 2019 are Albania, Greece, and the Czech Republic [31]. Of course, the role of imported cases is also important in new emergent diseases like Covid-19 where vaccination coverage is not present at all. An illustrative example of their importance is given by the data of Covid-19 in Iceland (as of June 3, 2020): 343 cases out of 1806 were reported to be infected abroad [32].
Here, we consider the occasional introduction of new cases into the population, at a very small rate ε, the rate of imported cases. From a mathematical point of view, these new cases will not destroy the deterministic oscillations obtained in Ref. [17] from a Hopf-bifurcation as long as ε is low enough and, at the same time, they prevent stochastic epidemic oscillations from extinction, as it happens when periodic solutions reach very low levels of prevalence (see Sec. IV C).
The main goal is, then, to analyze the SAUIS-ε model and to perform stochastic simulations of epidemics on networks to compare the regions of the parameter space where oscillations occur under both (deterministic and stochastic) approaches. Such a comparison will show that a significant part of the original region of the parameter space where deterministic oscillations occur is preserved when stochastic simulations are performed under the presence of a low rate of imported cases. Moreover, we will see that, as expected, a low level of awareness in the population at the moment of their introduction raises the chances of sparking sustained transmission and, hence, of generating new waves of infection. This could be, for instance, the situation in many countries once lockdown restrictions for Covid-19 have been lifted because most of their inhabitants are not immunized yet, as it happened in Singapore with a second wave of Covid-19 infections within its population of foreigner workers [33].

II. THE SAUIS-ε MODEL ON A GENERAL NETWORK
According to the SAUIS model introduced in Ref. [17], each node of a network of size N can be in one of the following four states: S (susceptible), A (aware), U (unwilling), and I (infected). Given any node n, let us denote the probabilities that n will have the respective states at time t by S n (t ), A n (t ), U n (t ), and I n (t ). Note that S n (t ) + A n (t ) + U n (t ) + I n (t ) = 1 for all time t 0, since each node has to be in one of the four states. Moreover, for the sake of simplicity, we assume along the paper that the alerting rates per link α 0 a , ν 0 a , α 0 i , the infection transmission rates per link β 0 , β 0 a , β 0 u , the decay rates δ a , δ u , and the recovery rate δ are the same for all the nodes. For instance, if at time t a node n is aware and one of its neighbors, m, is susceptible, then the probability that n successfully alerts m during the time interval (t, t + t ) is α a t + o( t ) provided that α a t < 1. Similarly, the probability that a noninfected node contracts the infection from abroad (imported case) during a time interval of length t is ε t + o( t ). Here it follows a summary of all transitions (or reactions) defining the SAUIS-ε model: Let us recall that the SAUIS model tries to account for the degradation of information quality among individuals that are aware of the epidemic situation. As well as in the standard SAIS models, the transition I + S → I + A represents the creation of a new aware individual that has acquired first-hand information about the epidemic by means of a direct contact. Also, the transition A + S → A + A creates new aware individuals that get indirect information from their acquaintances. Such new aware individuals have the same responsiveness as the information disseminators, which is not always the case. That is why the SAUIS model also includes the A + S → A + U transition, where U stands for unwilling to disseminate information. So, an unwilling individual has a lower level of awareness, in the sense that he or she does not try to convince other people about the risk, having in addition a weaker behavioral response.
The original SAUIS model in Ref. [17] contemplates the possibility that an infected individual, after recovering, may become aware with probability p, unwilling with probability q, as well as susceptible with probability 1 − p − q. Under these two additional transitions, periodic solutions are also possible (cf. Figs. 7 and 8 in Ref. [17]) but, since our ultimate goal is to provide evidence of robustness of the oscillatory regime in nondeterministic epidemics and to simplify the analysis, we will assume p = q = 0 along the paper. For the sake of simplicity, we have omitted them in the previous description of possible transitions.
Let us derive the approximate discrete-time equations for the evolution of these probabilities. An exact (but unfeasible) description would require the probability of the system being in any of the 4 N possible states. So, to derive approximate equations we will assume, as usual, that the joint probability for nodes n and m to be, respectively, in states X and Y is independent of the neighborhood's configuration of n and m. That is, it equals the product of both probabilities. This hypothesis will allow us to close the system without considering higher order terms for the joint probabilities (see Ref. [13] for a related discussion for the SAIS model).
Let t > 0 be small enough in such a way that, for every occurrence rate κ of a single event, the probability for this event to happen in the time interval (t, t + t ) is κ t + o( t ). For 1 n, m N, let a nm be the (n, m) element of the N × N adjacency matrix of the contact network, i.e., a nm = 1, if the nodes n and m are first neighbors, and a nm = 0 otherwise. With these ingredients, we can now write . The term multiplying A n (t ) corresponds to the event that the node n keeps being aware at time t + dt provided it was aware at time t, so it is 1 minus the sum of the probabilities of the three competing events that change the state A to another one: A → U with probability δ a dt, A → I with probability εdt, and the event that one or several infected neighbors of n succeed in performing the transition I + A → I + I (with probability β 0 a dt). For simplicity, the probability of this third event is computed as 1 minus the probability that none of such neighbors succeeds. Analogously, the term multiplying S n (t ) accounts for the probability that n is aware at time t + dt provided it was susceptible at time t. Now observe that, neglecting terms of order o(dt ), the expressions of the form m [1 − a nm κdtX m (t )] read as 1 − κdt m a nm X m (t ). This yields a dt m a nm I m (t ) + S n (t ) α 0 i dt m a nm I m (t ) + α 0 a dt m a nm A m (t ) , where in the second equality we have neglected again the terms of order o(dt ). Subtracting A n (t ) to both sides of the previous equation, dividing them by dt, and letting dt → 0, we obtain the differential equation governing the time evolution for A n (t ). Proceeding along the same lines for the other probabilities, we finally arrive at the following system of 3N ODEs: where the equation for S n (t ) is omitted because it is redundant (the sum of the nodal probabilities is always equal to 1). From the solution of Eq. (1) endowed with an initial condition, we can compute the expected number of aware, unwilling and infectious nodes at time t by summing the corresponding probabilities over the whole network, that is, Similar approaches to derive a system of equations for the probabilities for a node of being in one of several disease states have been previously introduced for the study of epidemics on networks and have received different names like, for instance, microscopic Monte Carlo approach in a discretetime setting [34], or N-intertwined model in a continuous-time setting [35]. An extension of the latter to multilayer networks is given in Ref. [36].

III. THE SAUIS-ε MODEL ON REGULAR RANDOM NETWORKS
To analyze Eq. (1) we start by the simplest case. So, let us consider the model over a random regular network (not necessarily fully connected) of degree k. As we will see, in this particular case the solutions of these 3N ODEs can be identified with the solutions of a much simpler system of three ODEs.
On this sort of networks, every node has the same vulnerability against the disease (in this setting, the degree is the only characteristic that distinguishes one node from another). So, it is reasonable to assume the same initial probabilities of being aware, A n (0) = a 0 , unwilling, U n (0) = u 0 , and infected, I n (0) = i 0 , for all nodes (what we call a uniform initial condition).
Now we focus on uniform solutions of Eq. (1), defined as those solutions A n (t ), U n (t ), I n (t ) that are independent from n. So, we can write A n (t ) = a(t ), U n (t ) = u(t ) and I n (t ) = i(t ) for all 1 n N. Then, the sums in Eq. (1) reduce to m a nm A m (t ) = k a(t ) and m a nm I m (t ) = k i(t ) because each node has the same degree k. So, the time evolution of these probabilities satisfies the following initial value problem (IVP): s + a + u + i = 1, endowed with the initial condition a(0) = a 0 , u(0) = u 0 , and i(0) = i 0 . Here all the alerting and infection transmission rates are per node and not per link, that is, Then, the local existence and uniqueness of solutions for both IVPs implies that a solution for Eq. (2) is also a solution (uniform by construction) for Eq. (1). That is, they are equivalent formulations of the epidemic on regular networks (see Lemma 3.1 in Ref. [13] for a proof of a similar result for the SAIS model).
Notice also that N A (t ) = a(t )N, N U (t ) = u(t )N, and N I (t ) = i(t )N. So, instead of thinking of nodal probabilities, we can consider the expected fraction of aware, unwilling and infected individuals as the macroscopic description (state variables) of our system which is more convenient to compare theoretical predictions with the outputs of the stochastic simulations of the epidemics.
As an incidental remark, it turns out that uniform equilibrium solutions for the complete system of 3N equations are only possible over regular networks. To see it, set A n = A, U n = U , I n = I for all 1 n N and denote the degree of any node n by k n . At the equilibrium, we have that 0 = dA n (t )/dt = k n (α 0 Considering a pair of different nodes n, m, from dA n (t )/dt = dA m (t )/dt we easily get that k n = k m . In consequence, all nodes must have the same degree.

A. Behavior of the equilibria
Similarly to the SAUIS model in Ref. [17], the tetrahedron := {(a, u, i) ∈ R 3 : 0 a + u + i 1} is positively invariant under the flow of Eq. (2). In fact, the vector field on the boundary of points strictly toward its interior for ε > 0. In consequence, there are no equilibria on the boundary of .
When ε = 0 the SAUIS model may have three different kinds of equilibria in : the trivial equilibrium P 1 = (0, 0, 0), the disease-free equilibrium Note that several distinct endemic equilibria may coexist. When ε > 0 in the SAUIS-ε model, the first two kinds of equilibria cannot exist and all possible equilibria are endemic. In fact, if ε is a small parameter the equilibria of SAUIS-ε can be interpreted as perturbations of the equilibria P 1 , P 2 and P 3 of the unperturbed SAUIS model. To study how these equilibria behave when ε is included in the model, we denote by f j (a, u, i; ε) the right-hand side of the jth equation in (2). Any equilibrium e * (ε) = (a * (ε), u * (ε), i * (ε)) of the system is implicitly given by the equations f j (e * (ε); ε) = 0, j = 1, 2, 3. We can derive implicitly the equations with respect to ε and get In the particular case of the equilibrium e * (0) = P 1 = (0, 0, 0) the previous system for ε = 0 can be solved and , In order that the trivial equilibrium P 1 of the SAUIS model stays inside the biologically feasible region when perturbed by ε, the three previous expressions must be positive. This occurs when β < δ and α a < δ a . In any other case, the SAUIS-ε model has no endemic equilibria bifurcating from P 1 for ε > 0 small. In fact, the eigenvalues of P 1 for the unperturbed SAUIS model are In consequence, P 1 bifurcates to an endemic equilibrium for the system SAUIS-ε when P 1 is hyperbolic stable in the SAUIS model. Moreover, by the hyperbolic property the equilibrium remains stable for ε > 0 small in the SAUIS-ε model. The basic reproduction numbers R 0 = β/δ and R a 0 := α a /δ a (see Ref. [17]) provide a clear interpretation of this fact: For the SAUIS-ε model with ε > 0 small enough, a stable equilibrium with all coordinates positive and small emerges when R 0 < 1 and R a 0 < 1, corresponding to a nonspreading, dying out SAUIS epidemic. In this case, the small equilibrium values of aware, unwilling and infected are essentially fed by the introduction of new infection cases at a rate ε rather than by the epidemic propagation itself.
When the perturbation is considered from the equilibrium Thus, the condition such that the perturbation of the diseasefree equilibria is inside the region for ε > 0 small is We point out that the expression on the left-hand side of the previous inequality is the same as the expression of the unique eigenvalue of P 2 for the unperturbed SAUIS model that may take positive values (see Eq. (13) in Ref. [17] and comments surrounding). The other two eigenvalues are either negative or have negative real part. In particular, this expression is negative if β < δ (since β > β a , β u ), meaning that di * dε | ε=0 is positive and, so, the SAUIS-ε system has an endemic equilibrium bifurcating from P 2 for ε > 0 small. When β > δ, the expression can be positive or negative depending on the other parameters. This change can be controlled by taking β a as a bifurcation parameter and so the bifurcation value is as shown in Ref. [17]. In the SAUIS model, the system shows a transcritical bifurcation as β a passes through the bifurcation value and the authors illustrate that this bifurcation may occur in two different directions. That is, by changing the stability of P 2 , from an stable equilibrium P 2 for β a < β c a a forward stable endemic equilibrium may bifurcate; or from an unstable equilibrium P 2 for β a > β c a a backward unstable endemic equilibrium may bifurcate. In both cases, the disease-free equilibrium is stable if β a < β c a and unstable otherwise. In consequence, an endemic equilibrium bifurcates from the disease-free equilibrium in the SAUIS-ε when P 2 is hyperbolic stable. As before in the case of P 1 , 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 toward its interior for ε > 0. This would be in contradiction with a hyperbolic unstable equilibrium entering from the boundary.
A similar treatment of the effect of immigration of infected individuals on the disease-free equilibrium of a general epidemic model (without awareness), in terms of the basic reproduction number, is given in Ref. [37].
Concerning the endemic equilibria P 3 , by means of the implicit function theorem, we know that the root e * (ε) of Eq. (2) will persist inside for ε > 0 small enough under the classical transversal condition as well as its stability. We refer to Ref. [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 particular, for ε = 0, a Hopf-bifurcation curve H 0 in the (β a , α i ) parameter space was obtained in Ref. [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 H 0 (see Fig. 1). For each value of ε, the regime of sustained oscillatory solutions of Eq. (2) lies inside the region of the parameter space limited by β a = 0 and the corresponding Hopf-bifurcation curve. Outside this region, solutions tend to a stable endemic equilibrium. This behavior is due to the fact that the real eigenvalue λ r of the Jacobian matrix J of Eq. (2) at the endemic equilibrium is always negative for any point of the considered region. Actually, λ r is smaller than the real part of the conjugate pair of complex eigenvalues λ ± that constitute, together with λ r , the spectrum σ of J. So, this means that the stability modulus of J, namely, max{Re(λ) | λ ∈ σ (J )}, is given by Re(λ ± ).
We recall that to compute the Hopf-bifurcation curve in the (β a , α i ) parameter space we need to find the solutions (β * a , α * i , a * , u * , i * ) of the system of equations given by the three equilibrium equations of Eq. (2) together with the condition that follows from Theorem 2.1 and Table 1 in Ref. [39] which guarantees that the Jacobian matrix J at the endemic equilibrium has a pair of pure imaginary eigenvalues. Precisely, this condition is where c 0 = − det(J ), c 1 is the sum of the principal minors of J, and c 2 = −trace(J ). From Fig. 1 we see that, for ε > 0, the region limited by H ε is entirely contained in that limited by H 0 . It follows that the more ε increases, the more the upper part of the original region of the oscillatory regime is reduced, whereas the lower part of this region remains almost the same for the three bifurcation curves. This fact suggests the existence of an abrupt transition from the nonoscillatory regime below the lower branch of the curves to the oscillatory one in the upper side of this branch. This transition is hardly perturbed by the external infection rates ε, provided that they are small enough.
In fact, the qualitative behavior of the solutions of Eq. (2) presents significant differences between parameter values near the upper and the lower boundary of the region determined by the Hopf-bifurcation curve. Although all solutions outside the confined zone are foci, a simple study of the eigenvalues of the Jacobian matrix at the endemic equilibrium shows that the transition from focus state (damped oscillations) to a periodic regime is faster through the lower boundary than through the upper one.
More precisely, it follows that the absolute value of the derivative of the stability modulus that controls the oscillatory motion is much larger for lower values of α i on the bifurcation curve (lower branch) than for higher ones (upper branch). In Fig. 2 we represent the stability modulus for β a = 0.1 in terms of α i . This fact becomes crucial to understand the difference of sensibility between the two boundaries in the stochastic Hopf-bifurcation diagrams appearing in next sections. A good example of this fact is the diagram presented in Fig. 9, where the simulations delimit the lower boundary of the Hopf curve almost identically as the theoretical curve. However, the upper boundary is hardly identified without taking into account the amplitude of the signal.

A. General simulation setup
As usual in the setting of continuous-time stochastic simulations, we use the well-known Gillespie algorithm on graphs [40]. All networks are randomly generated using the configuration model algorithm [41].

B. The special case of regular random networks
It is worth noticing that the classical Gillespie algorithm on graphs is highly time-consuming when executed over a network of about 10 4 nodes, in such a way that it is not feasible to construct a bifurcation diagram on two parameters, p 1 and p 2 , running 50 experiments for each pair (p 1 , p 2 ), when in addition the number of pairs is of the order of 10 3 . In the particular case of regular random networks, this serious drawback can be overcome by using what we will call the fast Gillespie algorithm (FGA in what follows). The FGA crucially relies on the following mean-field hypothesis (MFH): on a regular random network of big enough degree, the probability that a neighbor of a node has a given state can be approximated by the fraction of nodes on the entire network having that particular state. Let us see how the MFH can be used to speed up the Gillespie algorithm.
Assume that during an experiment over a given network a susceptible node n gets infected. In this case, the Gillespie algorithm updates the state of n (from S to I) and then explores all neighbors of n to update the number and type of links to be considered in the next time step. For instance, if a neighbor of n is aware, we lose a link of type A − S with associated weight α 0 a + ν 0 a , and we gain a link of type I − S with associated weight β 0 + α 0 i . This exploration of the neighbors of a node through an adjacency matrix (usually a pointer of pointers), that has to be done at every discrete time step, is one of the main computational loads of the classical Gillespie algorithm on graphs.
But assume now that the approximation given by the MFH assumption is good enough. Let N, N S , and N I 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 type I − S can be simply computed as kN I N S /(N − 1). Since the infection event I + S −→ I + I has rate β 0 = β/k, the total weight associated to all such events is then βN I N S / (N − 1). Analogously, the A particular event as, for instance, I + S −→ I + I is chosen with probability (N I N S β/(N − 1))/R. In this case, we just increase N I 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 and explore the neighbors of any particular node, simply because the MFH assumption allows us to work just with the absolute numbers of aware, infected, and unwilling nodes. Observe that FGA is independent of the degree k. In other words, k is not a parameter of the algorithm. The FGA performs statistically exact simulations of Markovian epidemic processes over regular random networks of high enough degree (that is, as long as the MFH applies). It is worth mentioning that what we have called FGA can be identified with the original version of the algorithm [40], which was aimed at the stochastic simulation of a fully mixed chemically reacting system. Let us see to which extent the outputs of the FGA and the Gillespie algorithm are essentially equivalent. Recall (Sec. III B) that for δ = 1, δ a = 0.01, δ u = 0.05, β = 3, β u = 0.5, α a = 0.01, ν a = 1, and ε = 10 −4 , a Hopf-bifurcation curve was obtained in the (β a , α i )-plane. In Fig. 3 we have shown the time evolution of the fraction of infected nodes according to the numerical integration of Eq. (2), together with the adaptive averaged outputs (see Sec. IV D) of the Gillespie algorithm over a regular random network of size N = 1000 and degree 50, and FGA with N = 1000. When α i = 0.04 and β a = 0.1 [ Fig. 3(a)], we are below the Hopf curve and so we have an stable endemic equilibrium. Observe that the three curves are essentially identical. For α i = 0.14 and β a = 0.1 [ Fig. 3(b)] we are inside the Hopf curve and we get a stable periodic orbit. In this case, the outputs from Gillespie and FGA seem qualitatively equivalent up to stochastic fluctuations.
Of course, we should give a precise meaning to the sentence seem qualitatively equivalent up to stochastic fluctuations. This is in fact the aim of Sec. IV D, where we give a detailed explanation about the statistical treatment of the simulation data to test the significance of the oscillatory regime. In Fig. 4 we show the complete Hopf-bifurcation diagram for the detection of the oscillatory regime in (β a , α i ) ∈ [0, 0.7] × [0, 1] after processing the output data obtained by both the classical Gillespie algorithm over a regular random network of N = 1000 nodes and degree 50 [ Fig. 4(a)] and the FGA with N = 1000 [ Fig. 4(b)]. We stress that this figure is intended only to compare both algorithms. In particular, we have chosen here N = 1000 since a higher order for N makes the computation of the Hopf-bifurcation diagrams under the Gillespie algorithm on networks highly costly 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 even for degrees as small as 50 (over 1000 nodes).

C. Impact of imported cases in the prevention of stochastic extinctions
Imported cases have an impact on the evolution of an epidemic if they happen when both disease prevalence and awareness level in the population are low. In their absence (ε = 0), an epidemic becomes eventually extinct because the number of infected nodes becomes extremely low when the level of awareness in the population is very high. In such a situation, it is well known that stochastic extinctions are extremely likely.
We can observe this fact in Fig. 5, which shows that all the region of the (β a , α i )-space where periodic solutions (interior of the Hopf-bifurcation curve) and weakly damped oscillatory solutions are predicted by Eq. (2) lies within the extinction zone (dark region). Note that, for a fixed value of α i , the greater β a is, the higher the disease prevalence at the endemic equilibrium and, hence, the lower the extinction probability of the epidemic. Conversely, for a fixed β a , the higher the alerting rate α i is, the lower the prevalence because aware individuals are more easily created and, hence, the higher the extinction probability is. This is the reason why the nonextinction region corresponds to the lower right part of this figure.  (β a , α i ). also, the moments where external infections are introduced. Such an oscillatory behavior is, in fact, predicted by Eq. (2). However, the extremely low prevalence attained in each cycle makes the stochastic extinction of the disease unavoidable, unless the occasional introduction of imported cases takes place when population awareness is low. The crucial role of such an introduction in maintaining the epidemic is revealed in Fig. 6(b), where the arrival of imported cases has been forced to cease (ε = 0) just at the beginning of the fourth flare-up. As expected, the fifth flare-up (dashed line) that would appear by keeping ε = 10 −4 now does not happen and the epidemic dies out. The total number of imported cases from t = 0 to t = 1200 is 119 in Fig. 6(a), whereas it is equal to 85 in Fig. 6 The low number of imported cases in the previous example, about 1 case every 10 infectious periods on average in a population of size 1000, shows that it is not necessary to have a high number of imported cases to prompt the occurrence of important flare-ups in populations whose individuals have a low level of awareness. This situation reminds, for instance, of what happened in New Zealand after the praised management of their first wave of COVID-19, where a reemergence of cases occurred after weeks with no community cases once social restrictions were lifted [42].

D. Statistical significance of the oscillatory regime
Our approach for the statistical detection of oscillations in the time series has two parts. The first one consists in the generation of three signals from the output of the stochastic simulations, one for each state variable of the model. As stated in Sec. IV A, for each combination of parameters we run K = 50 independent experiments. Each realisation consists on 2 12 points, of which we eliminate the first half seeking for stationarity of the time series, ending up with M = 2 11 points. In stochastic epidemic models, the usual way to construct a single signal from data is through the average. It is well known that ODE systems are a good approximation of the mean of realisations of stochastic processes in systems with a large number of components (nodes, molecules in reaction systems, etc.). However, although the standard average of trajectories works for endemic equilibria with a high prevalence [see Fig. 3(a)], it is not always a good option when a system exhibits dynamics where some of its components are present in low numbers [43].
In particular, the relationship between behavior and disease we are considering allows the existence of epidemics with periods of very low prevalence once individuals have become aware and adopt efficient protective measures against infection. Such periods are followed by relaxation in the adoption of preventive measures because of a decay in awareness which, in turn, is translated into new epidemic flare-ups. During these time intervals of low prevalence [see Fig. 6(a)], the stochastic trajectories of the epidemic show significant random fluctuations which lead to a phase shift in the time series with respect to the deterministic trajectory given by the solution of the ODE model. Consequently, when the average of trajectories is computed, the resulting signal has an important decrease of the amplitude of the oscillations and even a deformation of the periodic component.
To overcome this difficulty, we perform an adaptive mean of the signals. The idea is to locally align the signals before computing the average, which acts in favour of preserving periodic motion when the mean is performed. The first part of this method is classic in signal processing. To begin with, we perform a moving average computed over a sliding window of length 21 centered about each element in the time series. This acts as a low-pass filter to each time series attenuating the signal and omitting extreme frequencies. Then, taking the first time series, say X 1 , as a sample, we determine the time delay between X 1 and each of the other signals. To accomplish this, for each time series X k , k = 1, . . . , K, we compute the cross-correlation between X 1 and X k . The position τ k where the maximum is reached, that is corresponds to the position where the signals are best aligned. For each of the signals X k we define a new one,X k , by applying a circular shift of τ k positions to X k . In practice the previous cross-correlation is computed for small delays to ensure that the alignment is local. The adaptive mean we consider is given by The second part of the detection of oscillations consists in proving the statistical significance of the periodicity of the resulting data Y . To do so we define the periodogram of the zero-mean time seriesŶ := Y − Y by for n = 1, . . . , M 2 − 1, and with ω n = n M t and t being the time step of the time series. The periodogram is suitable to find hidden periodicities on the data. Indeed, data following a pure random process show a flat periodogram. However, if the time series presents a periodic motion, then the frequencies involved are shown as peaks of the periodogram. Fisher's 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 series contains a deterministic periodic component of unspecified frequency. The statistic taken into account to reject the null hypothesis is the maximum of the periodogram compared with the sum of all the values of the periodogram. That is, .
The previous maximum, as well as the summation, is taken for n = 1, . . . , M/2. This is so because the term P(ω 0 ) contains only information about the mean of the data. This term plays no role in the detection of oscillatory motion and, in fact, vanishes for zero-mean time series. The test is performed by computing the realised value g * of g from the data and then use the exact distribution of g given by Fisher in Ref. [44], where b is the largest integer less than 1/x. If the probability P(g > g * ) is less than α, then the null hypothesis can be rejected at level α. For further details on Fisher's test for hidden frequencies we refer to Refs. [45][46][47].
For the time series resulting from the parameters we are interested in, the direct application of the Fisher's test results positive for all configurations. This is so for two reasons. First, all motions associated to the solutions of Eq. (2) have an oscillatory component since they correspond to either a focus or a periodic orbit. Second, for the expression of the probability above, it easily follows that the larger the number of points M, the smaller the value of the probability. To be more precise in the determination of oscillatory motions, we reduce the points of the signals after the adaptive mean in a ratio 1 : 2 4 . That is, we end up with M = 2 7 points. This procedure enables to detect the strongest periodic motions in the parameter space, since weaker periodic signals will not pass Fisher's test with fewer points. Figure 7 illustrates it with two parameter configurations. On the left, periodograms with 2 11 points pass Fisher's test for both configurations. On the right, the periodograms of the same averaged signals with 2 7 points. In this case the configuration corresponding to the periodogram on the top side does not pass Fisher's test, whereas the one on the bottom side still passes the test. Notice that these parameter configurations correspond to the signals of top panels in Fig. 11. It is worth to mention that the frequencies avoided by the reduction of points are much larger than the frequencies of the data as can be seen in Fig. 7 and information about the periodicity of the signal is not lost.
For those parameters passing the Fisher's g-test with a p-value less than or equal to α = 0.01 we consider the corresponding peak frequency of the time series. The bifurcation diagram shows those parameters with an estimated frequency larger than the minimum observable frequency f min = 1/1500. A gradient of colors illustrates the amplitude approximation used in modeling of chemical reaction kinetics [48,49]. Such an approximation assumes that, in a system formed by different "chemical species" (nodes in different states in our context), the standard deviation of random fluctuations about the mean number of molecules of these chemical species scales as the square root of the size of the system (the number N of nodes in our network). So, dividing these mean numbers by N, it follows that, under this approximation, the standard deviation of the random fluctuations of the fractions of A, U , and I nodes is proportional to N −1/2 .

E. Oscillations in other network architectures
Another aspect of the robustness of the oscillations is their likelihood when other network topologies are considered. Are they still present when the epidemic spreads on more heterogeneous networks? If this is the case, then are they observed in the solutions of Eq. (1)? To answer these questions, we have used the configuration model algorithm to generate a Poisson network with mean degree 50 and an exponential network with mean degree 50 and minimum degree 25. The first degree distribution corresponds to the well-known Erdös-Rényi random graphs, whereas exponential degree distributions have been observed in empirical contact networks [50] and have a much higher variability.
Unfortunately, in this case it is not possible to reduce Eq. (1) to a simpler system, as it was done in Sec. III B for regular random networks. So, the algebraic approach used there for computing Hopf-bifurcation curves is no longer feasible. Instead, using β a and α i as tuning parameters (with increments of size 0.002 or even smaller when approaching the turning point of the curve), we have obtained an approximation to these curves by numerically integrating Eq. (1) for N = 1000 (i.e., the full system of 3000 equations) until a long enough time T 2 . Then, for each pair (β a , α i ), we computed the coefficient of variation (CV ) of the fraction of aware, unwilling and infected nodes from the solution A i (t ), U i (t ), and I i (t ) of Eq. (1) for the last 1000 units of time (T 2 − T 1 = 1000). Precisely, for the fraction Finally, an approximate Hopf-bifurcation curve is obtained by using the criterion that, for a given pair (β 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 periodic when we integrate the system up to a given time can become nonperiodic when longer times are considered. This is particularly relevant when trying to delimit the upper branch of the bifurcation curve. Here the transition from periodic solutions to very weakly damped solutions is hardly noticed because of its flatness, which was already observed when computing the algebraic curve for the regular case (see Fig. 2). To calibrate our criterion with respect to the integration time, we compared the algebraic Hopf-bifurcation curve for Eq. (2) with those obtained from Eq. (1) for a regular random network of degree 50. Notice that Eq. (2) and, hence, the corresponding algebraic Hopf-bifurcation curve, do not depend on the degree, while Eq. (1) does depend on it. Figure 8 shows that, as expected, increasing the integration time from T 2 = 10 000 to T 2 = 20 000 leads to a better approximation to the whole algebraic curve (solid line). This figure also shows that, for T 2 = 20 000, 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 Eq. (1) using an integration time T 2 = 20 000 and a threshold value of CV equal to 0.1% for the fraction of the three types of nonsusceptible nodes.

F. Simulation results and discussion
The values of the parameters used in the simulations are taken from Ref. [17] (except for the rate of imported cases ε) and reflect what we consider it is a natural scenario, although they are not intended to model any particular disease. First, aware (and unwilling) individuals are affected by lower transmission rates because of the adoption of preventive measures. Second, the mean infectious 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 without willingness to convince others about the risk of infection. 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 unwillingness period (1/δ u = 20).
The procedure of the adaptive mean of the signals from the stochastic simulations gives an excellent reconstruction of the oscillatory behavior of single trajectories. In Fig. 11, we show the comparison of a solution of Eq. (2), the adaptive mean of the 50 realisations of an epidemic on a regular random network of size N = 10 000 generated by means of the FGA, and the trajectory used as the reference for the alignment of the signals. Note that this alignment preserves the periodicity of the original signals (trajectories) but, at the same time, it reduces their amplitudes. With this respect, it is worth recalling that the amplitudes shown in the Hopf diagrams do not correspond to stochastic simulations of single epidemics, but to the adaptive means of 50 stochastic trajectories for each pair (β a , α i ) of parameter values (see Sec. IV D for details).
Hopf-bifurcation diagrams show that the abrupt transition between a nonoscillatory regime and an oscillatory one given by the lower branch of the Hopf-bifurcation curve is quite accurately captured in all the settings we have considered. The robustness of this transition is, in fact, already observed when this curve is computed for different rates of imported cases (see Fig. 1).
The best agreement is achieved for regular random networks of size N = 10 000. Figure 9 shows that the region with periodic solutions with amplitudes larger than 0.05 falls neatly within the limits predicted by the Hopf-bifurcation curve. Comparing this figure with the panels in Fig. 4 obtained for N = 1000, one can see that the agreement of the observed region of periodic solutions with high amplitudes (yellow and orange squares) with 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 enough for weakly damped oscillations to disappear.
A similar agreement with the predictions as the one in Fig. 4 is also observed in networks with different topologies but with the same number of nodes, N = 1000 (see Fig. 10). Moreover, we can observe that the differences between the Hopf-bifurcation curves when the degree distributions have different variability are remarkable. For regular random networks and Poisson networks, these differences are small, i.e., both types of networks lead to similar (but not equal) Hopf-bifurcation curves. However, the resulting curve for an exponential network is clearly different, with a turning point at a higher value of α i and lower values of β a [see Fig. 10(b)]. Interestingly, the region of the (β a , α i ) parameter space corresponding to periodic solutions with large amplitudes (yellow squares and most of the orange ones) clearly falls inside the region limited by the Hopf-bifurcation curve.
All in all, stochastic simulations confirm the predictions of the ODE model showing that the proposed mechanism for prompting awareness is able to generate recurrent epidemic cycles on different network architectures, including the one with an exponential degree distribution which has been claimed to describe contact patterns in real populations. Interestingly, in all cases the predicted region of the parameter space for the oscillatory regime satisfy that β a < β u , a natural condition according to the higher level of alertness assumed for aware individuals with respect to the unwilling ones. Moreover, these regions always contain those averaged stochastic trajectories with the highest amplitudes.

V. CONCLUSIONS
The interplay between human behavior and epidemic spread has been considered in many papers during the last 15 years [2,51]. Most of them have as the main goal to elucidate the dependence of the basic reproduction number on the different behavioral responses as, for instance, social distancing, rewiring of connections, awareness, etc. [10,11,13,14,23,25]. Few of them address the existence of oscillating epidemics [20,21,[27][28][29] and, as far as we know, only in Ref. [17] sustained oscillations arise uniquely from the interaction between awareness dissemination and epidemic spread, i.e., without the need of any recruitment of new (susceptible) individuals.
In this paper we have challenged the existence of oscillations predicted by the SAUIS model introduced in Ref. [17] by means of stochastic simulations on different network topologies. The model has been formulated on networks and considers the existence of a very small inflow of imported cases. These cases are essential to keep the epidemic going on because disease prevalence attains very low levels when there is a high degree of awareness in an oscillating epidemic. This is not a problem at all in a deterministic framework, but it leads to an unavoidable stochastic epidemic extinction in all parameter combinations where sustained oscillations are present.
Our simulations show that the presence of a small number of imported cases allows, when the disease prevalence is low and the number of aware and unwilling individuals decreases, an oscillatory behavior of the stochastic epidemics which closely resembles the one predicted by the deterministic counterpart of the model. In particular, the existence of an abrupt transition from a stationary regime where oscillations are strongly damped to an oscillatory one with amplitudes of more than 10% in the number of infected nodes predicted by the ODE model is clearly observed in the simulations over different types of networks. Moreover, when the size of the network increases (N = 10 000), periodic epidemics with amplitudes larger than 5% carried out on regular random networks fall within the oscillatory regime in the parameter space predicted by the ODE model (see Fig. 9). So, the robustness of the predictions about the existence of oscillations due to pure behavioral changes has been established.
An interesting example of oscillating epidemics is given by the evolution of the incidence rate of STDs during the last decades. The reemergence of STDs like syphilis and gonor-rhea occurring since the mid-1990s [52] has been associated with a decrease in awareness after the introduction of the antiretroviral therapy for HIV and, indeed, it appeared after an incidence decline in the 1980s. This decline coincided with the emergence of the global AIDS pandemic and has been attributed to preventive behavioral changes in response to HIV campaigns during that time [53]. However, it has been also claimed that the drop in incidence before 1984 occurred too early to be ascribed to such induced behavioral changes and may be part of a long term cyclic trend of this type of diseases (although the nature of this periodic behavior remains unexplained) [54].