key: cord-0687186-eko2zkul
authors: Ochab, Magdalena; Manfredi, Piero; Puszynski, Krzysztof; d’Onofrio, Alberto
title: Multiple epidemic waves as the outcome of stochastic SIR epidemics with behavioral responses: a hybrid modeling approach
date: 2022-03-15
journal: Nonlinear Dyn
DOI: 10.1007/s11071-022-07317-6
sha: 02add0f6fa7d4f0a1a5cf2705510d0de1130e512
doc_id: 687186
cord_uid: eko2zkul

In the behavioral epidemiology (BE) of infectious diseases, little theoretical effort seems to have been devoted to understand the possible effects of individuals’ behavioral responses during an epidemic outbreak in small populations. To fill this gap, here we first build general, behavior implicit, SIR epidemic models including behavioral responses and set them within the framework of nonlinear feedback control theory. Second, we provide a thorough investigation of the effects of different types of agents’ behavioral responses for the dynamics of hybrid stochastic SIR outbreak models. In the proposed model, the stochastic discrete dynamics of infection spread is combined with a continuous model describing the agents’ delayed behavioral response. The delay reflects the memory mechanisms with which individuals enact protective behavior based on past data on the epidemic course. This results in a stochastic hybrid system with time-varying transition probabilities. To simulate such system, we extend Gillespie’s classic stochastic simulation algorithm by developing analytical formulas valid for our classes of models. The algorithm is used to simulate a number of stochastic behavioral models and to classify the effects of different types of agents’ behavioral responses. In particular this work focuses on the effects of the structure of the response function and of the form of the temporal distribution of such response. Among the various results, we stress the appearance of multiple, stochastic epidemic waves triggered by the delayed behavioral response of individuals.

tic discrete dynamics of infection spread is combined with a continuous model describing the agents' delayed behavioral response. The delay reflects the memory mechanisms with which individuals enact protective behavior based on past data on the epidemic course. This results in a stochastic hybrid system with timevarying transition probabilities. To simulate such system, we extend Gillespie's classic stochastic simulation algorithm by developing analytical formulas valid for our classes of models. The algorithm is used to simulate a number of stochastic behavioral models and to classify the effects of different types of agents' behavioral responses. In particular this work focuses on the effects of the structure of the response function and of the form of the temporal distribution of such response. Among the various results, we stress the appearance of multiple, stochastic epidemic waves triggered by the delayed behavioral response of individuals.

Keywords Multiple epidemic waves · Social distancing · Hybrid systems · Stochastic epidemic models · Human behavior · Memory effects · Gillespie algorithm

The principles of mathematical epidemiology (ME) were established more than one century ago with the milestone work by Sir Ronald Ross in 1916 and by Kermack and McKendrick [1, 2] .

The Kermack and McKendrick's lied in the description of infection transmission by means of the mass action law of Chemistry and Statistical Mechanics. In their formulation, the individuals' social contact patterns relevant for transmission are represented as collisions between particles of a perfect gas, while transmission is modeled as a chemical reaction occurring with a certain probability upon a random encounter between two individuals of different types, namely a susceptible subject and an infective one [3] :

Susceptible + Infectious β → Infectious + Infectious.

As a consequence, in their approaches, the social contact rate per individual and the transmission rate β per social contact were taken essentially as natural constants valid for appropriate combinations of spatial settings, human activities, cultural habits, institutions, etc. However, it was acknowledged the possible dependence of these constants on the seasonality of social phenomena, e.g., the school calendar, or on climatic phenomena [4] .

In the last four decades, mathematical models of infectious diseases have become key supporting tools for public health decisions [5] [6] [7] . This development was made possible by the increasing awareness of the need for better data and statistical tools allowing a finer and finer description of the transmission process. Further dimensions were incorporated in the models: individual-level characteristics, as age and sex; mesolevel fundamental structures, such as the community composition by household type, by space, by social activity, etc. [5] .

The impulse impressed by the fear for avian influenza and the 2009 H1N1 flu pandemic induced the development of highly sophisticated mathematical and computational models, which better integrate models with data [6, 8] , and which were used to describe and predict the worldwide dynamics of influenza pandemics. Of course, the current COVID-19 pandemic, and related mitigation measures, are bringing further astonishing momentum to the discipline, with an endless list of contributions (more than 9000 modeling papers and PhD theses on COVID-19 are reported by Google Scholar at February 2, 2022).

However, in essentially all public health models, the individuals' social, transmission and vaccination behavior remain unaffected regardless of the trends in the risks of acquiring the infection they might perceive from available information.

Clearly, this is an unrealistic abstraction, increasingly less applicable in contemporary scenarios. This became even more true during the current COVID-19 pandemic, which represented a huge open-sky laboratory of humans' behaviors [9] . Indeed, individuals are frantic information seekers and therefore hardly unaffected by the state of the disease and related information.

Attempts to remove this abstraction have led, in the last 20 years, to the birth of a new branch of ME that was termed the behavioral epidemiology (BE) of infectious diseases [7, [10] [11] [12] .

The pioneering work in BE was due to Capasso and Serio [13, 14] , who first assumed that the transmission intensity might also depend on the level of infection spread by modeling the contact rate β as a decreasing function of infection prevalence. This was the first instance of an epidemiological model incorporating individuals' spontaneous social-distancing behavior by a prevalence-dependent contact rate [7, [10] [11] [12] .

A key postulate of BE is that infection spread is not anymore an independent process, but it is instead the outcome of the interplay of an infection-spread layer with a number of other layers, first of all a behavior layer, in turn modulated by an information layer. In simple words, human behavior relevant for infection spread may be critically affected by the information about the spread and severity of the considered infection. And any change in human behavior will possibly affect infection trends that might eventually feedback on behavior itself. This means that the spread of information is a critical component of infection transmission modeling and therefore requires careful modeling per se. All this led some of us to develop new classes of epidemiological models by using new variables, that we termed information indexes, in order to include not only individual perceptions about the current state of infection, but also the memory of past spread [7, 15, 16] and eventually also spatio-temporal memories. In particular, in [17] the Capasso-Serio model was extended by considering a deterministic SIR model for an endemic infection disease with a general information-dependent transmission rate allowing to also include past available information on infection spread.

In this article, we depart from the previously cited works [13, 17] along two main directions.

The first one deals with the interpretation side: we will cast the aforementioned behavior implicit models including the effects of delayed information [17] ) in the framework of the theory of nonlinear feedback control systems [18] [19] [20] [21] . The underlying idea is that the individuals' behavioral responses considered in the BE literature, whereby, e.g., agents reduce their daily number of social contacts depending on the information they have on disease severity. This is nothing other than a (nonlinear) feedback control. In particular, unlike many feedback systems studied in mathematical biology [19, 22, 23] , the one considered here is a genuine feedback. Indeed, it is the outcome of voluntarythough uncoordinated-actions of agents.

The second main departure point is that most BE studies are based on deterministic models and as such they are suitable to describe epidemic or endemic scenarios for large populations only.

More in detail, we aim at studying the impact of human behavior during the spread of an epidemic outbreak in small-medium size populations, by taking into account that human decisions might also depend on past information. This has a number of implications. First, the involved nonlinear feedback control systems are stochastic and therefore depart from the classical asymptotic deterministic description based on the qualitative theory of differential equations [18] [19] [20] [21] . Moreover, as we are dealing here with a system having a finite life span, it is necessary to assess the impact of the control action on a mainly transient phenomenon.

In particular, we depart from classical stochastic models of mathematical epidemiology, generically belonging (as the models of chemical kinetics, of which they can be considered as a subcase [7] ) to the class of nonlinear birth and death Markov processes [24, 25] . The latter processes can be simulated by means of the well-known Gillespie algorithm [26, 27] . However, the presence of the information memory index (denoted as M(t) in this paper) makes, from the theoretical viewpoint, the resulting model a hybrid stochastic system [28] [29] [30] , where a birth-death model of epidemic spread is coupled with a piecewise deterministic model for the information index M(t). The main implications are the following: (i) the state variables of the stochastic epidemic sub-system impact on the-otherwise deterministic-dynamics of the memory; (ii) in its turn, the model of the memory impacts on the birth-death epidemic component by making the key occurrence probabilities, namely the transmission rate β, heterogeneous in time. From the computational viewpoint, the latter issue requires to modify Gillespie's algorithm following the lines of [28] [29] [30] [31] .

Based on this background, we then use selected simulation experiments to analyze and classify the effects of different individuals' delayed behavioral responses during the course of a stochastic SIR-type epidemic in a fully susceptible population. In particular, the agents' responses to the epidemic threat are modulated by: (i) the structure of the behavioral response function, (ii) the form and amplitude of the response function temporal kernels. Consistently, for each hypothesis on the behavioral response function, the resulting system is simulated to adequately reconstruct a number of key features of the stochastic epidemic, namely the simulated distributions of (i) the final attack rate, (ii) the epidemic extinction time, (iii) the epidemic maximum prevalence, etc.

Finally, we explicitly stress that our work did not intend to refer to COVID-19, which is an infection with a more complex dynamical structure than the one considered here [32] [33] [34] . Instead, our work had a quite general purpose, namely to develop a general stochastic behavioral epidemiology framework including delayed responses by individuals. However, we will insertonly where appropriate-observations that could be related to (or of interest for) investigations on the current pandemics.

The manuscript is structured as follows. Section 2 (i) presents the behavioral deterministic models and the related information indices that motivated this study, (ii) it plugs them in the formalism of nonlinear control theory, and (iii) it presents the models main features. Section 3 introduces the stochastic counterparts of the given models with special focus on the simulation issues that arise when lagged behavioral components are considered, thereby leading to hybrid stochastic systems. Section 4 explains the experiments. The subsequent sections illustrate the results for different forms of the information index. Concluding remarks follow.

In this section we present the deterministic counterparts of the stochastic behavioral models used in this paper. We depart from two basic, behavior implicit, models with a prevalence-dependent response. Then, we discuss the structure of the information index M and present some main instances of delaying kernels tuning the behavioral response, namely, the exponen-tially fading kernel and a generalization that we term the acquisition fading kernel. Then we set the proposed systems in the framework of nonlinear control theory. Finally, we discuss the main qualitative features of the proposed class of models with special focus on the interplay between the strength of the behavioral response and the related time delay.

2. 

where S = S(t), I = I (t) represent the number of susceptible and infective individuals at time t, γ > 0 is the recovery rate, and β(I ) is the prevalence-dependent transmission rate which was assumed to obey

where β 0 represents the transmission rate in the epidemic early phase where no individuals' responses are in place yet. The key hypothesis is that the larger the infective prevalence, the stronger will be the agents' response to the epidemic threat by correspondingly reducing their transmission-or contact-rates, that is by expanding their social distancing. Among the main properties of model (1)-(2) we mention that: (i) an epidemic outbreak will occur only if β(0) (S(0)/N ) > γ , i.e., if the basic reproduction number (BRN) relevant for this model exceeds one; (ii) for t → ∞ it is I (t) → 0, that is, due to the deterministic continuous nature of the model the outbreak only ends for infinitely times. Moreover, the total number of infected subjects obeys the biologically meaningful relationship

. In substantive terms, the working of the model is rather simple: the inclusion of a simple agents' response increasing with prevalence has a protective role on the community. Compared to the baseline Kermack-McKendrick model, this response can mitigate the epidemic in the event of an outbreak and can-if very intense-even prevent the outbreak occurrence. This simplicity will obviously disappear in presence of more articulated behavioral response.

In [17] , the notion of prevalence-dependent social distancing first formulated by Capasso and Serio was extended in a fairly general manner by taking the transmission rate β as a decreasing function of a general information index denoted by M(t), i.e., β = β(M), β (M) < 0. The concept of information index was first introduced in [15] , to summarize the entire amount of-current and past-available information about infection spread that agents can use to inform their response to the infection threat. In [17] the idea of social distancing was cast within a basic behavior implicit information-dependent SIR model for an endemic-rather than epidemic-infection (therefore also including the vital dynamics of the population). The model read as follows:

where μ represent both the birth and death rates and kept equal to ensure population stationarity over time.

The corresponding epidemic variant is obtained by just setting μ = 0. A model using the M index in a behavior-explicit setting was proposed by [35] within their outbreak model with normal vs altered behavior.

As argued in [17] , the information index is a flexible summary of the infection state and can therefore take a wide range of forms. A general form adequate for our present purposes is the following prevalence-dependent form

where g(I ) is an increasing function of the prevalence, and K (q) is a delaying kernel (or 'memory function'), that is a probability density function that weighting the levels of past prevalence. For the sake of the simplicity, in our subsequent analyses we will use g(I ) = k I . In [17] it was shown that in the endemic case (i.e., for μ > 0), model (3)-(4)-(5) has two equilibria: (i) the disease-free equilibrium point D F E = (N , 0), which is Locally Asymptotically Stable (LAS) and Globally Asymptotically Stable (GAS) provided that the model BRN is lower or at most equal to one, where BRN = β(0) μ + ν and a unique endemic equilibrium

where (S e , I e ) is the unique solution of the following system

. The stability of the endemic state depends on the specific model of information index (5) . For example, if the kernel K obeys K (q) = δ(q) where δ(.) is the Dirac function then M(t) = I (t) (i.e., the model is unlagged), and E E is GAS (see [17, 36] and Appendix A).

Even on an intuitive ground, the information index M(t) can be read as a feedback control variable (traditionally denoted as u(t) in nonlinear control theory), since it represents a decision variable that-at the individual level-agents may tune to reduce their hazard of infection. By setting β 0 = β(0) and defining (5) in the classical form of nonlinear control theory [18] x

where

In particular, the case where the information index is assumed to be simply given by the infection prevalence, M = I = x 2 , has a nice control-theory interpretation: it is the deviation of the state x 2 (t) (i.e., I ) with respect to the ideal reference case of absence of disease, i.e., x re f

In this section, we state our main assumptions on the transmission rate and the delaying kernel, and provide some relevant control-theoretic interpretations.

As regards the functional form of the transmission rate resulting from the above definition (2.3), i.e., β = β 0 (1 − ψ(M)), this could allow to design the control by applying techniques of differential geometry [18, 20, 21] . However, this approach would not keep under control the key constraint on function ψ(M), namely:

Therefore, in our numerical simulations we will use the following phenomenological family of functions for the information-dependent transmission rate β(M):

where p is an integer such that p ≥ 1, and M 50 is such that β(M 50 ) = 0.5β(0) Note that, in the particular case where the information index only includes current information (that is, M(t) = I (t)) and for p = 1, the corresponding force of infection (FOI) h(I ) = β(I )I /N (the instantaneous hazard of infection faced by a susceptible individual per time unit) is increasing, whereas for p ≥ 2 the FOI is nonmonotone with a maximum at I = M 50 .

Note moreover that, the fact that the DFE is LAS (GAS, actually) only if BRN < 1 independently of the memory kernel , i.e., only if in absence of the behavioral effect the disease cannot remain endemic, implies that one cannot design a ψ(M) such that the DFE is stabilized.

As regards the memory kernel, we will consider three main cases:

• The memoryless case K (q) = δ(q). In this case the information index M only includes current information.

• The exponentially fading memory kernel (EFK):

This represents the often realistic case where past information receives a decreasing weight. In this case, the integro-differential system can be reduced to ordinary differential equation (ODEs) by adding the further ODE:

• The acquisition-fading memory kernel (AFK):

This represents another possibly realistic scenario where no information (about the infection state) is available at current time; then, it is gradually acquired up to a maximum before starting fading out exponentially, as in the EFK case. Even in this case integro-differential model (3)-(4)- (5) can be reduced to ODEs by the additional equations

Note that if the two sub-processes of information acquisition and memory fading occur at the same relative speed a 1 = a 2 = a, then our 2-dimensional acquisition-fading kernel reduces (via the limit a 2 → a 1 ) to the so-called Erlang(2, a) = ate −at kernel [37] .

It is immediate to verify that:

• In the case of the exponentially fading kernel, the memory index/control variable M(t) can be read as the output of a low-pass filter with characteristic cutoff frequency a that is applied to the 'signal' g(I (t)). Under our work assumption g(I ) = k I the transfer function H (λ) between the input I and the output M reads as follows

where λ ∈ C and M(λ), M(λ) and K (λ) are the Laplace transforms of, respectively, M(t), I (t) and K (q); • In the case of the acquisition-fading kernel, the memory index M(t) can be read as the output of the application of two low-pass filters in series (with characteristic frequencies a 1 and a 2 ) that is applied to the 'signal' g(I (t)). In particular, in the case g(I ) = k I , the transfer function between the input I and the output M reads as follows

Therefore, in our problems, the processes of fading and (if it is the case) acquisition of the memory are analogous to the error with memory connected to the measure of a feedback signal by means of a mechanic or electronic device [19] . This is a typical scenario in control theory, where the controller often acquires the information on the variables to be controlled by means of measure devices. The latter can have their own internal dynamics and, moreover, are also aimed at reducing extrinsic noises that are characterized by large frequencies.

Adopting again the control-theory formalism, the final structure of the considered family of models (6) is as follows:

• In case of EFK:

• In case of AFK:

We discuss here the main stability properties of the endemic state of model (6) under the different hypotheses considered on the form of the delaying kernel. First, in the case of EFK, it was proven in [17] that the endemic equilibrium is LAS. Moreover, in [38] by means of an appropriate Lyapunov function it was shown that under the additional condition (Mβ(M)) > 0 the endemic equilibrium is globally stable. Additionally, in [17] , the case of an Erlang 2 -type memory was considered. As already pointed out, this memory represents a special case of the acquisition-fading kernel where the two sub-processes of acquisition and fading occur at the same rate a 1 = a 2 = a. The corresponding transfer function reads as follows

In this case, the EE can be either LAS or can be destabilized by a Hopf bifurcation (see [17] ).

Let us now consider model (6) for a general delaying kernel K (τ ). Denoting as K (λ) the corresponding Laplace transform and linearizing at EE one gets:

Under an AFK, K (λ) is given by formula (11) and it can easily be rewritten as follows:

This finally yields the following characteristic polynomial:

Since all the coefficients are positive, applying the Routh-Hurwitz criterion we get the following local stability conditions of the EE:

In particular, if

then EE becomes unstable through Hopf a bifurcation at the locus

As well known, the Hopf theorem only provides local information in the parametric space. In other words, it only states that in a parametric neighborhood of the bifurcation point the solution of the system is a periodic solution of small amplitude. Therefore, in our case this will hold for values of ( p, σ ) in the region H ( p, σ ) < 0 that are sufficiently close to points where H ( p, σ ) = 0. In order to draw more general information about points in the interior of the region H ( p, σ ) < 0 one can take advantage of the Yakubovitch-Efimov-Fradkov theorem [39] .

To apply the theorem, note preliminarily that the system is bounded and has two unstable equilibrium points, and the DFE has as stable manifold the segment l = (x 1 , 0, 0, 0) with x 1 ∈ [0, N ] that corresponds to the trivial case of absence of disease in the target population. The biological and geometrical nature of the stable manifold of DFE implies that there are not heteroclinic orbits connecting EE to DFE.

As a consequence, from the Yakubovitch-Efimov-Fradkov theorem (see corollary 1 of [39] ) if H ( p, σ ) < 0 then system undergoes self-sustained oscillations, not necessarily periodic, called Yakubovitch oscillations or Y-oscillations [40] . We report below two useful remarks.

Remark The above results can easily be extended to delays in the form of n ≥ 3 first-order low-pass filters in series.

Remark Note that, in our case (and in the generalization) the Y-oscillations are fully induced by the behavior, i.e., by the control. In control language, this is the consequence of the fact that the control is based on a 'measure' coming from a double filtering of the prevalence, i.e., from the indirect nature of the control. Note however that the choice of the weight K (q) is primarily based on behavioral arguments.

As regards the analytical form of condition (12) , this is cumbersome and gives no particular hints, so that also (13) becomes of scarce practical relevance. Therefore, we numerically computed the Jacobian and its 4 associated eigenvalues at the endemic equilibrium. This allowed us to plot (see Fig. 1 ) the stability and instability regions in function of the two key parameters: (i) T 1 = 1/a 1 which is the average acquisition time of information; and T 2 = 1/a 2 which is the average fading time. In particular, Fig. 1 reports the contour plot of the following function:

In the rest of the present work we will mostly focus on the purely epidemic case (μ = 0). This case was not previously analyzed with the exception of the case M = I , i.e., the Capasso-Serio model. For the general model with a generic informationdependent M index it is easy to see that it exists a disease-free equilibrium x 2 = 0 and that x 1 (t) tends to an equilibrium value that depends on the function ψ(M(t)) and on the dynamics of M(t), which is in turn determined by the delay kernel. A fact not previously pointed out in the literature is that under appropriate behavioral responses multiple behavior-induced epidemic waves can arise even in the course of a single epidemic. We will numerically illustrate these deterministic results in "Appendix."

In this section we introduce the stochastic counterparts of the previous deterministic models with social distancing and discuss the simulation issues that arise when the behavioral component is involved. We start from the simpler case where no memory effects are considered, and then extend to the case of hybrid systems arising under memory effects.

As is well known, the deterministic models of infectious diseases can be considered as a good approximation of the true underlying stochastic systems only if the population size is large enough so that stochastic oscillations around the mean pattern become negligible. In the absence of behavioral effects, if the population is small then the system is modeled by a birth-death Markov process [24] .

In the absence of memory effects, i.e., when considering the stochastic version of the Capasso-Serio model, the resulting stochastic model remains a birthdeath stochastic processes [24, 26, 27, 41] .

In this case, the model is simply the outcome of two transition processes and related events, namely: (i) transmission: a susceptible subject becomes infectious due to adequate social contacts with infectious subjects. The probability of a transmission event during the infinitesimal interval (t, t + dt) is:

However, unlike the basic stochastic SIR model, the probability in formula (14) is not bilinear. (ii) removal: an infectious individuals are removed from the infectious compartment due to recovery and immunity acquisition. The probability of a removal event during (t, t + dt) is:

The stochastic simulation of birth and death Markov processes is elegantly done by the Gillespie algorithm [26, 27, 41] which in our case works as follows. Let us suppose that the nth event occurred at time t n so that the system after the event has state values

Therefore, denoting the time of the next event as

then τ is such that

where r n is exponentially distributed. This yields:

The type of the next event, which will either be the contagion of a susceptible or the removal of an infectious subject), is determined by randomization through a Bernoulli experiment based on the two complementary probabilities

In the presence of a general memory effects the behavior of the system becomes more complicated. Indeed, the stochastic evolution of the discrete state variables (S, I ) has to be coupled with an integral equation describing the dynamics of M(t), which is a continuous one during each inter-events intervals (t n , t n+1 ).

In particular, under our hypotheses on the memory kernels, the dynamics of M is given by a system of ordinary differential equations depending on I n . As a consequence, the system is no longer a classical birth and death process, but it belongs to the class of the so-called hybrid models [28] where stochastic automata are coupled with deterministic processes by means of time-varying propensities. Note however that in our case the coupling is mutual. Indeed, a stochastic birth-death model (for infection spread) is coupled with an ODE (or a system of ODEs) model for the memory, yielding for M a continuous piecewise deterministic model that depends on the stochastic epidemic model. Indeed, in between two stochastic events, the dynamics of M(t) is fully deterministic, i.e., M(t) is a piecewise deterministic process.

Overall, the above-described model can be formalized similarly to the class of stochastic hybrid automata defined in [31] .

Operatively, the infinitesimal probability of a contagion event in (t, t + dt), given by:

is time-dependent. This requires to resort to a version of the Gillespie algorithm with time-varying propensities, as in [28] . Namely, the time to the (n+1)-th event is determined by the following equation:

In the general case, the inequality

As a consequence, from (18) and provided that I n > 0 (the case I n = 0 is trivial because I = 0 is an adsorbing state for the system) we can infer that

where τ L = r n γ I n + β 0 I n N S n τ R = r n γ I n which is an useful bound for the numerical simulations.

Remark The bounds τ L and τ R are not constant: they change at each step of the simulation algorithm.

To compute τ it is necessary to know the model linking the continuous state variable M(t) to the discrete state variables (S, I ) .

Once τ has been determined, the type of the next event is chosen (as in the classical Gillespie algorithm) by the following probabilities:

Then one has to numerically solve the following equation γ I n τ + I n N S n ψ(τ ) = r n (19) in the interval

by some numerical algorithms. We then have two main cases. The first one occurs when function ψ(τ ) is available in closed form. In the case where β(M(t)) can be analytically integrated in (t n , t n+1 ), then Eq. (19) can be solved by means of well-known standard iterative algorithms, such as the Newton-Raphson and the bisection methods (see "Appendix"). For the sake of notational simplicity we define f (τ ) = γ I n τ + I n N S n ψ(τ ) − r n thus to determine τ we must solve

Note also that: f (τ L ) < 0 and f (τ R ) > 0. Instead, if ψ(τ ) is not analytically known, one has to numerically approximate the integral that appears in the time-varying Gillespie algorithm. For this case, we propose the algorithm described by the following steps:

• Split the interval [0, τ R ] in L >> 1 points, and define

Note that:

• Approximate ψ(θ k+1 ) by some standard strategy, such as:

+β(M(t n + θ k+1 ))) .

Note that this implies that 

is the integer part of x, e.g., [3.99] = 3.

• Suppose that the step q ≥ k L had been reached and that

Then go to the next step θ q+1 and compute f (θ q+1 ) by means of formula (20) . The following three cases are possible: (i) | f (θ q+1 )| < T ol, where T ol is a sufficiently small tolerance, so you can stop; (ii) f (θ q+1 ) < −T ol, thus one must continue; (iii) f (θ q+1 ) > T ol, which implies that then the solution τ is such that τ ∈ (θ q , θ q+1 )

Since h is sufficiently small, one can use the approximation

The most sensitive point in the previous algorithm is therefore that of appropriately choosing V , which must be sufficiently large. In our simulations, we have chosen V in a way such that the time step h is always kept smaller than 0.1 days. Moreover, one must remember that since τ R changes at each step of the Gillespie algorithm, the same holds for V Once determined τ the type of the next event is chosen by the following probability: ρ rec = γ I n γ I n + β(M(t n + τ )) I n N S n ρ cont = 1 −ρ rec 3.3 The case of the exponentially fading memory

In the cases of EFK and of AFK kernels, the coupling between the discrete and the continuous parts of the hybrid model is relatively simpler than in the general case. Indeed, as we have seen, in our examples the model of M(t) reduces, respectively, to one or two linear ordinary differential equations.

In this section and in the following one, we will show how β(M(t)) can be analytically defined in all event intervals (t n , t n+1 ) for both the EFK and the AFK. For the exponentially fading memory kernel, we will also provide the analytical form of ψ(τ ), which allows to significantly speed up the computations.

In the case of an exponentially fading memory, one has that in each interval (t n , t n+1 ) between two consecutive events, M(t) is given by the following linear differential equation

since we suppose that the epidemics starts at t = 0. The solution of ODE problem (21)- (22) is:

where

The function ψ(τ ) is given by the integral

As for the above integral note the following: (i) if σ = 1 then we have that the upper integration limit is smaller than the lower integration limit (due to a > 0), but since y = R + exp(−aθ) ⇒ y > R, thus the function to be integrated is negative and the integral is positive; (ii) if σ = −1, the upper integration limit is greater than the lower integration limit and since y = R−exp(−aθ) ⇒ y < R, the integrand function turns out to be positive, so that the integral will be positive as well (and, in principle, analytically solvable). Further setting y = K z one gets:

where r = R/K and ρ = σ/K . Namely, for the case p = 1 one has:

or equivalently (the following version reduces numerical errors in the case exp(−aτ ) << 1)

For the case p = 2 one has:

For p ≥ 3, the expression of ψ(τ ) becomes untractable.

Of course, once τ has been computed, one must set M n+1 = g(I n ) + (M n − g(I n ))Exp(−aτ )

In the case of the acquisition-fading kernel, one has that in each interval (t n , t n+1 ) the solution (Z (t), M(t)) to the pair of addition differential equations can be computed as follows. First, Z (t) is defined as the solution of the following differential equation

The solution of the above ODE problem is: 

implying:

Note that in this case the integral defining ψ(τ ) cannot be analytically solved even in the case p = 1.

Remark In this case, once τ has been numerically computed, one must set

Simulations of scenarios where the memory kernel is the AFK are shown in appendix, for the sake of the readability. In our work we used MATLAB 2019b for model simulation and the analysis of results. Given our focus on multiple-stochastic-epidemic waves, we briefly present the approach for the localization of (multiple) epidemic peaks. The purpose of the algorithm is that of allowing to separate true 'process' peaks, i.e., peaks induced by the behavioral response, from spurious oscillations due to the stochastic nature of the problem. The localization of the first peak was done by finding the maximum number of infectious individuals, because the first peak always resulted the higher one. The latter (purely due to susceptible depletion) is a somewhat straightforward consequence of the sim-plicity of the model. The localization of the second peak (and subsequent ones) was more complicated due to the noise of stochastic simulations. Firstly, we filtered the time courses and analyzed only cases with 3 and more points. Noise was filtered by the Gaussianweighted moving average filter (we applied MATLAB function: smoothdata) with window size equal to 100 points. Then we accepted as true peaks, maxima arising in the smoothed data course according to the following criteria (i) minimum height of peak equal to 100 infectious, (ii) minimum time distance between subsequent peaks of 5 days (MATLAB function: findpeaks). As matter of fact, since the appearance of multiple peaks was much neater in the case BRN = 15 (compared to BRN = 2), later on we will mainly report our results on multiple peaks for this case. 

In this section we report the main results of a range of stochastic simulation experiments, with special focus on the effects of different hypotheses on the shape of the social distancing response function β(M) on the system behavior.

In all our experiments we simulated epidemic outbreaks in a fully susceptible population, under the following assumptions on input parameters: (i) the population size is set to N = 10,000; (ii) the infection average duration is set to one week, implying γ = 1/7days −1 ; (iii) M 50% is allowed to take three possible values: M 50% ∈ {10, 50, 200}. As regards the basic reproduction number, we considered two widely different cases:

(1) BRN = 15 corresponding to a highly transmissible infection, such as measles; (2) BRN = 2, corresponding to a moderately transmissible infection as it may be the case, e.g., of seasonal influenza. The corresponding baseline transmission rates (in the absence of any behavioral response) are β 0 = 15/7days (−1) , and β 0 = 2/7days (−1) , respectively. In the deterministic SIR epidemic model in a wholly susceptible population, the considered figures of the BRN would yield a total proportion of people eventually infected, i.e., a final attack rate of the outbreak-well above 99percent for BRN = 2. As for the main outputs of stochastic simulation we primarily looked at the simulated probability density functions (PDF) of (i) the final attack rate, (ii) the epidemic extinction time, (iii) the maximum prevalence, and to the temporal realizations (or realized paths) of the prevalence function. In particular, in order to adequately reproduce the relevant probability density functions, the model simulation is replicated N sim = 100,000 times.

In this subsection, we considered the memoryless case with p = 1, in order to characterize the role of parameter M 50% . The stationary PDFs of the final attack rate and of the extinction time are shown in the upper and lower panels of Fig. 2 Concerning the final attack rates (upper panels of 2), we note that for both BRN = 15 and BRN = 2, the simulated PDFs are always bimodal. The mode for low values of the infected subjects reflects stochastic extinction occurring soonafter infection introduction, when the behavioral effects are still negligible. Obviously stochastic extinction is rather unlikely for BRN = 15, whereas it occurs with a much higher probability percentages for the 'classical SIR'; instead for BRN = 2 this lower mode is very large.

We may note that even for M 50% = 200 for both BRN = 15 and BRN = 2 the effect is remarkable. Qualitatively, both for BRN = 15 and for BRN = 2, the mean value of the total number of people that acquired the infection decreases if M 50 decreases, but the variance of the distributions increases. For BRN = 2, at M 50 = 10 the distribution is characterized by a large initial peak followed by exponential-like decrease.

As for the extinction time (lower panels of 2), the simulated PDFs are also bimodal, corresponding to "initial" and "final" (i.e., after the outbreak occurred) extinction. Interestingly, for BRN = 15 the mean values of the extinction time are increasing as M 50 decreases (i.e., when outbreaks occur, the stronger the behavioral response the longer is the epidemic), and also the variances become larger. Moreover, at M 50 = 10 the PDF is more markedly bimodal. For BRN = 2 and M 50% = 10 the PDF again shows an initial peak followed by exponential-like decrease.

Note that, for BRN = 15 the location of PDFs of the extinction time for M 50% = 50 and especially for M 50% = 10 is so rightward that the chosen approximation of neglecting vital dynamics could be questioned. The same holds for BRN = 2 and M 50% = 50. This suggests that behavioral responses could significantly delay the epidemic extinction time and, not paradoxically, sustain the endemicity of the infection. Both these effects have been observed during the current COVID-19 pandemic.

For the sake of completeness, the right panels of Figures 2 and 3 report analogous results for the case BRN = 2. Note in particular the much larger stochasticity arising compared to BRN = 15, a well-known fact.

The PDFs of the maximum reached prevalence and of the time of occurrence of such maximum are reported in the upper and the lower panels of Fig. 3 , respec- 

The effects of the shape of the behavioral response function β(M), in absence of delays, is summarized by parameter p. Of course, even the threshold M 50% plays a fundamental role. Figure 4 reports a sample of the simulated time series of the infectious prevalence with M 50% = 200, for different values of p (we keep BRN = 15). This shows three main effects of increasing p. The first is related to the fact that, for p > 1, the function β(I ; p) decreases slower with I , thereby increasing extinction time. This effect is better illustrated by the corresponding PDFs in Fig. 5 . The second and the third are more interesting: (i) the maximum epidemic size decreases with p and (ii) the disease prevalence remains at large values for a larger time w.r.t. the case p = 1. Particularly interesting is the case p = 100, whose behavior can be better understood by resorting to the following approximation of function β(I )

Initially, since for I < M 50% it holds β(I ) ≈ β(0), the epidemic grows free up to the level I = M 50% . Sooner or later I switches to the value M 50% where β = 0.5β(0). However, sooner or later the next contagion event will occur, so that prevalence I reaches the value M 50% +1 and the transmission rate abruptly falls to β = 0. This means that the contagion stops and prevalence will remain constant at the level I = M 50% + 1 until the next removal event occurs, so that the system now switches back and forth around M 50% . In other words, I (t) will experience for a long time small stochastic oscillations around I = M 50% until when through a series of removal events disease extinction is reached.

By decreasing M 50% the minimal extinction time remarkably increases. This is better shown in Fig. 5 . Figures 6 and 7 report the results for the case BRN = 2. As for the total number of infected people, the corresponding simulated PDF gradually switches rightward as p increases (see supplementary figures).

In this section we investigate the effects of delayed social distancing response function β(M) on the epidemic course. In particular, we will study the EFK case by considering two possible values for the parameter a tuning the delayed response: a = 0.1/day and We assessed the effects of a delayed behavioral response in both the stochastic and in the deterministic case (see "Appendix"), for BRN = 15 (Fig. 8 ) and BRN = 2 (Fig. 9 ) and for distinct values of p ( p = 1, 2, 10). Both deterministic and stochastic simulations show that multiple epidemic waves can actually occur, as a consequence of the behavioral response of individuals, when lagged information is used in the response process. The different cases reported in the figures clarify mechanisms and determinants of these waves. Notably, these waves cannot occur in the absence of the delay. In particular, the oscillatory pattern will be more marked the larger the amplitude of the time delay (i.e., the lower a), the sudden the behavioral response (i.e., the larger p), and the larger the response threshold M 50% .

In what follows, we reports the main statistics of our simulative result where multiple waves occur. To do so, in the computation of the simulated PDFs we did not include those realizations where no second peak was clearly detectable. In the present experiments we keep M 50% = 50, whereas a ∈ {0.05, 0.1}, p = 1, 2, 10, BRN ∈ {2, 15}. Figure 10 shows the PDF of infectious prevalence at the first epidemic peak. We note that of course in all cases there is a peak at low prevalence levels, corresponding to a rapid extinction of the epidemic. A second peak occurs for larger values of prevalence. We note that for BRN = 2 (bottom part of the figure) the PDFs are located at low values of prevalence, they are overlapping and the average values are increasing with p. On the contrary, for BRN = 15 (top part of the Figure) , the PDFs are located at large values of the prevalence, they are not overlapping and the average values are decreasing with p. Obviously, the larger the delay in the behavioral response, the larger the expected magnitude at the first peak. Figure 11 reports the PDF of the infectious prevalence at the second epidemic peak. As detailed before, this is done only for the case of a large value of the BRN (BRN = 15), in order to magnify occurrence and detectability of the second wave. For a = 0.1 the prevalence at the second time series peak decreases with p, whereas for BRN = 15 and a = 0.05 the his-tograms overlap, but their location appears somewhat nonmonotone with p.

As for the simulated PDF of the time of occurrence of the first epidemic peak (Fig. 12 ) , we note that (disregarding the initial peak due to stochastic extinction) for BRN = 15, the PDFs are well-separated with an average peak-time decreasing in p (from 4 to 8 days), while for BRN = 2, the PDFs largely overlap especially for longer behavioral response delays. In particular, the effects of the response delay are essentially negligible on the timing of the first peak.

Instead, the corresponding PDFs of the time at the second epidemic peak (still drawn for BRN = 15 only) reveal ( Fig. 13 ) a marked effect of the time delay in the behavioral response.

We report in Table 3 the main summary statistics of the PDFs shown in the previous figures.

For the sake of completeness, we report also the PDFs of the time to epidemic extinction (Fig. 14) . Still disregarding initial stochastic extinction, it can be noted that no relevant role was played by the delay in the behavioral response. However, the dependency on p and the BRN remains more complicated.

As for the epidemic final attack rate, the corresponding simulated PDFs (Fig. 15 ), suggest-as usual-a major role of the shape of the behavioral response (tuned by p) for BRN = 2. For BRN = 15 a counterintuitive role of the delay emerges at high values of p: for p = 10 and a = 0.05, the PDF of the final attack rate has two nontrivial modes. One of these is the traditional one at high levels of the FAR. The second one occurs for levels of the FAR around value 20%. This second mode reflects the interesting fact by which a too delayed behavioral response cannot prevent the occurrence of the outbreak but, if the response is sufficiently intense, it can halt it and bringing it to extinction.

Finally, the PDF of the number of nontrivial peaks (Fig. 16) in the more interesting case (BRN = 15) , reveals the number of (detectable) peaks increases with the delay of the behavioral response (that is, for a = 0.05 compared to a = 0.1).

For BRN = 2, Fig. 17 compares, for a single realization of the stochastic process, the simulated time series of the information index M(t) and that of the infectious prevalence I (t). The figure well shows the analogy between the concept of an exponentially fading memory and the concept of 'low-pass filter' used in system control theory. Calling these quantities as the (original) 'signal' (I (t)) and the 'filtered signal' (M(t)), it happens that the signal is noisy and the filtered signal is smooth, for a = 0.01 and for (a, p) = (0.1, 1), whereas in all the other cases also the filtered signal is noisy, especially for a = 0.05. This is a natural consequence of the fact that the 'low-pass frequency' a is too close to the typical frequency of the stochastic events of I (t). Case BRN = 15 is shown in Appendix.

This article investigated the effects of agents' behavioral responses on the dynamics of stochastic SIR models for epidemic outbreaks. The agents' behavioral response was represented by means of an informationdependent transmission rate specified as a function of an appropriate information index, as first proposed in [15, 17] . In particular, agents were assumed to respond either to the current or past prevalence of infection, where the latter was specified according either to an exponentially fading memory or to an acquisitionfading memory. This resulted in a family of hybrid stochastic models. In order to numerically simulate specific models belonging to the family, we needed a suitable extension of classical Gillespie algorithm. In particular, we developed analytical formulas valid for two classes of models.

Our analysis provides a thorough classification of the possible outcomes of behaviorally modulated stochastic SIR epidemics depending, besides basic reproduction numbers, on: (i) the form and strength of the behavioral response, (ii) the information time lag with which this behavioral response is enacted by individuals.

This offers a range of theoretical results. For example, even in the absence of delays in the behavioral response (i.e., when the agents' response is only based on current prevalence), the behavioral response may substantially increase the variance of the extinction time, thereby favoring infection persistence (especially for highly transmissible infections). This behaviortriggered persistence might link with other factors, such as demographics or the appearance of new variants of the pathogen, favoring endemicity, as might sadly be the case with COVID-19. These phenomena becomes richer when the behavioral response is also modulated by past information. In this case, we showed that multiple epidemic waves can result from delayed agents' responses already in the underlying deterministic models and remain also in the stochastic formulation.

In general the topic of recurrent waves, especially in relation to endemic infectious diseases (e.g., for measles), has been a major topic of mathematical epidemiology over decades [6] . Instead, the issue of multiple waves during an epidemic outbreak has become popular during the pandemic preparedness in the early 2000s. This started from retrospective studies on the 1918-1919 Spanish flu pandemic and developed after the H1N1 2009 pandemic. A variety of explanations were proposed, ranging from exogenous changes in transmissibility [42] , to policy intervention [43] , to the presence of multiple strains [44] , up to behavioral changes [45] . Recently, elegant analytical work has stressed the role of inter-regional commuting in inducing multiple (behavior-unrelated) epidemic waves [46] . The multiple waves observed during the COVID-19 pandemic in industrialized countries were primarily due to the activation and relaxation of social distancing measures during the first pandemic year and to the onset of new variants of concern during 2021 [47], but also behavioral explanations were proposed [48] .

In relation to the previous debate, the main innovation of this work lies in the fact that it represents the first general theoretical assessment of the role of delayed behavioral responses in promoting multiple waves in the most basic stochastic setting. This in turn required to consider an innovative formal setting, namely that of Stochastic Hybrid Automata, whose continuous component mirrors (in our case) a distributed delay/memory.

Of course, the proposed model is minimalistic: many refinements can be considered. First, the model represents spontaneous individual behaviors in the simplest manner, i.e., implicitly by the information index approach [15, 17] . As such, the model might also represent a situation where a government uses the same information index to implement a policy of forced behavior change in individuals as, e.g., a lockdown. Modeling behavior in an explicit manner [7, 15, 16] can refine such scenarios.

Further improvements might include: (i) further epidemiological detail, e.g., the presence of a latency period; (ii) the effects of exogenous seasonal fluctuations in transmission; (iii) the combination of individuals' behavior with governmental responses; (iv) the spatial dimension modeled as continuous or by a discrete metapopulation approach ( [46] ): it is well-known that the metapopulation network topology can relevantly impact on disease spread [7, 49] ; (v) the presence of multiple pathogen strains; (vi) the possibility that, besides the transmission rate, the behavioral response affects also the recovery rate. For COVID-19 this surely occurred, e.g., in the propensity to testing and selfisolating by pauci-symptomatic individuals; (vii) the modeling of the propagation of information through its articulated interpersonal and virtual networks [49] [50] [51] . δ(x) without imposing technical conditions. In paper [36] , it was shown that the endemic equilibrium is globally stable provided that F(I ) = Iβ(I ) is either increasing or having at most an inflexion point. Here we remove this requirement.

Since model (1) ) an appropriately chosen topologically equivalent dynamical system (TEDS) obtained by multiplying the velocity field by a Dulac multiplier [52] . In our case, we employed the following multiplier: 1/(Iβ(I )). The resulting TEDS reads as follows

where we set for the sake of notation simplicity: N = 1 , A(I ) = 1/β(I ), and c = γ + μ. The rescaled time τ is linked to the original time t by τ = τ 0 I (t)β(I (t))dt.

Note that at the endemic equilibrium it is which implies that at equilibrium the Jacobi matrix J is such that The output of the deterministic model for p = 1 is reported in Fig. 18 , where the effects of three different values of parameter M 50% are shown.

B.2 The effects of the shape of the response function in the deterministic case For sake of comparison with our stochastic simulation, in Fig. 19 the behavior of the deterministic model in the memoryless case is reported for two values of p.

Here we illustrate the effects of a delayed behavioral response in the deterministic case, for BRN = 15 (Fig. 20) and BRN = 2 (Fig. 21) and for distinct values of p ( p = 1, 2, 10). In this section, we will illustrate the results of the simulations in case of acquisition fading kernel. Since the acquisition phase is short, we set a 1 = 1, whose characteristic time is T 1 = 1. On the contrary, the fading phase is much slower so we considered two values: a 2 = 0.1, i.e., T 2 = 10 >> T 1 , and a 2 = 0.05, i.e., T 2 = 20 >> T 1 . Similarly, multiple waves are also observed for BRN = 2 (not shown). Patterns of multiple waves are also observed in the stochastic model, for both BRN = 15 (Fig. 24 ) and BRN = 2 (Fig. 25 ). Figure 26 shows the prevalence at the first peak. For BRN = 15 it is interesting to note as the location of the PDFs depends nonmonotonically on p.

As far as the prevalence at the second peak, Fig. 27 shows that for a 2 = 0.1 the PDF for p = 1 overlaps with the PDF for p = 2, which for a 2 = 0.05 moves rightward and overlaps partially with the PDF for p = 10.

As far as the time of the first peak, Fig. 28 shows that the PDFs seem mixtures between an unimodal PDF and an exponential PDF (faster decaying for R B N = 15 , slower for BRN = 2).

As far as the time of at the second peak, Fig. 29 shows that for BRN = 15 the average decreases with p and also the variance.

The distribution of the time at epidemic extinction time is of interest (see Fig. 30 ) . Indeed for (BRN, a 2 ) = (15, 0.1) the location of the PDFs is increasing with p, but the PDF for p = 10 has also an appreciable earlier mode. On the contrary, for (BRN, a 2 ) = (15, 0.05) all the three PDF have an earlier peak, which in the case p = 10 is the largest one. For BRN = 2 the PDFs overlaps are mixtured between unimodal and fastly decaying exponentials.

Finally, even the distribution of the final attack rates is of some interest (Fig. 31) . Indeed, for (BRN, a 2 ) = (15, 0.1) the PDF are increasing with p, but for p = 10 the PDF has a smaller but appreciable earlier mode. For (BRN, a 2 ) = (15, 0.05), all the three PDFs have an intermediate mode, and for p = 10 this mode is the largest one.

Finally, Fig. 32 shows as p is a key factor influencing the average and shape of the histogram of main peaks. 51. Wang, Z., Xia, C.: Co-evolution spreading of multiple information and epidemics on two-layered networks under the influence of mass media. Nonlinear Dyn. 102(4), 3039-3052 (2020) 52. Wiggins, S., Wiggins, S., Golubitsky, M.: Introduction to Applied Nonlinear Dynamical Systems and Chaos. Springer, Berlin (1990) Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

An application of the theory of probabilities to the study of a priori pathometry-Part I

A contribution to the mathematical theory of epidemics

Mathematical Structures of Epidemic Systems

Seasonality in epidemic models: a literature review

Infectious Diseases of Humans: Dynamics and Control

Modeling Infectious Diseases in Humans and Animals

Statistical physics of vaccination

Comparing large-scale computational approaches to epidemic modeling: agent-based versus structured metapopulation models

Non-pharmaceutical interventions during the covid-19 pandemic: a review

Modelling the influence of human behaviour on the spread of infectious diseases: a review

Behavioral epidemiology of infectious diseases: an overview

Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases

A generalization of the Kermack-Mckendrick deterministic epidemic model

Global solution for a diffusive nonlinear deterministic epidemic model

Vaccinating behaviour, information, and the dynamics of sir vaccine preventable diseases. Theor

Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases

Information-related changes in contact patterns may trigger oscillations in the endemic prevalence of infectious diseases

Nonlinear Control Systems: Analysis and Design

Feedback for physicists: a tutorial essay on control

Nonlinear Systems Analysis

Analysis and Design of Nonlinear Control Systems

Feedback Control in Systems Biology

Mathematical Physiology

The Elements of Stochastic Processes with Applications to the Natural Sciences

Stochastic Methods

A general method for numerically simulating the stochastic time evolution of coupled chemical reactions

Exact stochastic simulation of coupled chemical reactions

Adaptive simulation of hybrid stochastic and deterministic models for biochemical systems

Tumour suppression by immune system through stochastic oscillations

The interplay of intrinsic and extrinsic bounded noises in biomolecular networks

Distributed delays in a hybrid model of tumor-immune system interplay

Spread and dynamics of the covid-19 epidemic in Italy: effects of emergency containment measures

Effects of information-induced behavioural changes during the COVID-19 lockdowns: the case of Italy

The geography of covid-19 spread in Italy and implications for the relaxation of confinement measures

The effect of risk perception on the 2009 h1n1 pandemic influenza dynamics

Epidemic models with nonlinear infection forces

Time Lags in Biological Models

Global stability of infectious disease models with contact rate as a function of prevalence index

Oscillatority of nonlinear systems with static feedback

Conditions of existence of oscillations for hybrid systems

Markov Processes: An Introduction for Physical Scientists

Coinfection can trigger multiple pandemic waves

Qualitative analysis of the level of cross-protection between epidemic waves of the 1918-1919 influenza pandemic

Increased transmissibility explains the third wave of infection by the 2009 h1n1 pandemic virus in England

The effect of public health measures on the 1918 influenza pandemic in US cities

Exactly solvable sir models, their extensions and their application to sensitive pandemic forecasting

Awarenessdriven behavior changes can shift the shape of epidemics away from peaks and toward plateaus, shoulders, and oscillations

An sir model with infection delay and propagation vector in complex networks

A new propagation model coupling the offline and online social networks

The authors warmly thank three anonymous referees of this Journal whose comments allowed to greatly improve the exposition of the manuscript. Usual disclaimers apply.

Since part of this work has been done during the stay of Dr Ochab at the international Prevention Research Institute (IPRI), a private firm, the authorization to get the simulated time series obtained during her stay at IPRI must be asked to IPRI. Code availability Since part of this work has been done during the stay of Dr Ochab at IPRI, a private firm, the authorization to get the codes prepared at IPRI must be asked to IPRI. However, this work contains a detailed description of the algorithms allowing all readers to replicate their implementation.

The authors declare that they have no conflict of interest.

In this appendix we show the global stability of the endemic equilibrium in the no-memory case K (x) =