Hysteresis loop of nonperiodic outbreaks of recurrent epidemics.

Most of the studies on epidemics so far have focused on the growing phase, such as how an epidemic spreads and what are the conditions for an epidemic to break out in a variety of cases. However, we discover from real data that on a large scale, the spread of an epidemic is in fact a recurrent event with distinctive growing and recovering phases, i.e., a hysteresis loop. We show here that the hysteresis loop can be reproduced in epidemic models provided that the infectious rate is adiabatically increased or decreased before the system reaches its stationary state. Two ways to the hysteresis loop are revealed, which is helpful in understanding the mechanics of infections in real evolution. Moreover, a theoretical analysis is presented to explain the mechanism of the hysteresis loop.

INTRODUCTION

The spreading of epidemics is currently one of the hottest topics in the field of complex networks, and a great deal of significant progress has been achieved so far, including the infinitesimal threshold [1-6], the reaction-diffusion model [7-10], flow-driven epidemics [11-15], objective spreading [16,17], temporal and/or multilayered networks [18-27], and other aspects [28-36]; see Refs. [37-39] for details. A common point in all these contributions is that their infectious rate remains constant during an evolutionary process, i.e., will not change for the whole evolutionary process from an initial infected seed to the final stationary state. It is well known at present that there is a critical , although it may be different for different cases. In the thermodynamic limit, the infected fraction will be zero when , and it will become nonzero at . After the critical point , the infected fraction will gradually increase with a further increase of . Thus, a key question for all these cases is how to determine . To figure out the answer, a bifurcation diagram of is usually calculated by using the same or similar initial conditions for each , i.e., one or a few initial infected seeds. To be specific, for each in the bifurcation diagram we let the system evolve until it attains its stationary state, and then we record the stabilized infectious density [37]. Then, we turn to a new and perform the same process. Thus, the obtained bifurcation diagram of describes how the stabilized infectious density depends on the parameter . In this sense, we are concerned with how an epidemic breaks out, i.e., the growing phase where changes from to .

However, we noticed from real data that on a large scale, the spread of an epidemic is in fact a recurrent event with both growing and recovering phases [40-42]. That is, there is a recovering phase with gradually decreasing , and a growing phase with gradually increasing . This recovering phase is distinctive from the recovery process at the individual level with fixed . For example, in the traditional SIS model, will remain unchanged when the infected individuals recover to susceptible status. Additionally, in the SIR model, will also remain unchanged in the evolutionary process in which the infected fraction will first gradually increase to its maximum and then decrease until zero. A primary purpose of studying epidemic spreading is to understand it and then to control or eradicate it. For the same reason, it is also necessary to understand the recovery phase, but little attention has been paid to this subject so far.

On the other hand, it was pointed out that the infectious rate is in fact dependent on time and is seriously influenced by seasonal variation [40-45]. To reflect this influence, Ref. [40] assumed that depends sinusoidally on time , and Ref. [45] assumed that depends on the infected density. Based on these assumptions, an interesting finding is that it is possible for the transition from to to be a first-order phase transition with hysteresis [44,45], in contrast to the second-order transition observed in all the cases of constant  [37-39].

We recheck several real data here [46,47], and we pay attention to both the features of growing and recovering phases and the features of a transition with hysteresis. We find that (i) the growing and recovering phases do not overlap but take different paths, forming a hysteresis loop (see Fig. 2 as an example), and (ii) there is no jumping behavior in the hysteresis loop, indicating that the hysteresis is not a first-order but rather a second-order transition. To reveal the underlying mechanism of this new hysteresis loop, we present an epidemic model that incorporates two characteristic features of real data. The first one is the influence of seasonal variation. A seasonal epidemic spreading generally lasts for a few months, i.e., beginning from a healthy state, then gradually developing into an outbreak, and finally recovering to the healthy state again. In this process, the weather and environment will undergo an approximately periodic change, implying that will increase accordingly from zero to its maximum and then decrease to zero again. The second feature is the succession of status during the evolutionary process of epidemic spreading, indicating that the system's parameters are adiabatically changed.

FIG. 2.

A typical hysteresis loop in real data on epidemics. (a) Amplification of the data around the peak between 1938 and 1939 in Fig. 1(b), where the dashed line indicates the location of the rescaled zero point. (b) The hysteresis loop of (a) in the rescaled framework, where the squares and circles denote the growing and recovering phases, respectively.

Through the long history of epidemic study, much data have been collected and many epidemic models have been presented to explain the data [48]. These models have been adequate in explaining traditional epidemics. However, for modern epidemics, much evidence has shown that they spread on complex networks. For example, the data on gonorrhoea told us that the epidemic spreading paths are of network topology [49]. It has been widely confirmed in recent years that the topology of a network plays an important role in epidemic spreading [37]. Accordingly, the developed models and theories can be used in practical applications. For example, the original seeds of SARS can be successfully tracked by a small fraction of epidemic data [50]. In this sense, we will take here a typical complex network as the infrastructure for an epidemic to spread.

Using our model, we show that the hysteresis loop can be reproduced if is adiabatically increased or decreased before the system reaches its stabilized state. This way of changing in a transient process is more similar to the process of epidemic spreading in reality, in contrast to the previous ways in which is allowed to update only when the system reaches its stationary state, and the initial infected seeds for the updated will be reset [37]. Thus the critical role of the practical mechanics of infections in evolution is revealed. Moreover, we present a theoretical analysis to explain the underlying mechanism.

The paper is organized as follows. In Sec. II, we first analyze the real data and then introduce our model. In Sec. III, we perform numerical simulations to regenerate the phenomenon of the hysteresis loop revealed in Sec. II. In Sec. IV, we present a mean field to explain the numerical results. In Sec. V, we present conclusions and discussions. Finally, in the Appendixes, we present two more accurate theoretical approaches to explain the numerical results. One is the theoretical analysis based on the microscopic Markov-chain approach, and the other is the theoretical analysis based on the effective degree Markov-chain approach.

MODEL DESCRIPTION

By following the approach from data to model, we recheck here the recurrent epidemic data of Hong Kong, New York, and Baltimore [46,47], and we replot them in Figs. 1(a), 1(b), and 1(c), respectively. From these we see that each time series consists of a small amplitude background with occasionally large bursts, implying that an epidemic outbreak is unlikely to occur unless some necessary conditions, such as the match of seasonal weather, humidity, and sunlight, are satisfied. We also see from Figs. 1(a)–1(c) that each epidemic outbreak lasts a relatively long time, i.e., from one to a few months, indicating the mark of a season. On the other hand, we have to point out that each time series in Fig. 1 is not from one or a few initial seeds in the same initial stage, rather each time series is a sum from different initial seeds at different initial stages. For example, the data in Fig. 1(a) were collected by the Hong Kong Department of Health from a sentinel surveillance system with 50 general practitioners [46], i.e., a summation of 50 subtime series. In this sense, we conclude that the data reflect in fact the match between epidemic spreading and seasonal environment. A better match implies a better environment for epidemic spreading. That is, we may assume that the degree of match is closely related to the value of . Therefore, a burst (or an epidemic outbreak) in Fig. 1 can be modeled by a change of from to .

FIG. 1.

Real data of epidemics. Part (a) represents the weekly consultation rates of influenza-like illness (per 1000 consultations) in Hong Kong from the general practitioners (GPs) [46], where the value of is from 0 to 150. Parts (b) and (c) represent the time series of reported cases of infective measles in New York and Baltimore, respectively, where the variable is from 0 to in (b) and from 0 to in (c) [47].

Moreover, we notice from Fig. 1 that each burst consists of both growing and recovering phases, separated by a maximum of an infected fraction, i.e., a peak. For the former, the epidemic increases gradually until it reaches its maximum, while for the latter, the epidemic decreases gradually until it recovers to a healthy society again. To see this more clearly, we replot one typical burst of Fig. 1 in Fig. 2(a). With the help of the dashed line, we see that the growing and recovering phases in Fig. 2(a) are not symmetric. Let be the time of the peak in Fig. 2(a) and . Figure 2(b) shows the dependence of on , where the squares and circles denote the growing and recovering phases, respectively. It is easy to see that there is a hysteresis loop in Fig. 2(b).

If we ignore the aspect of how the data are collected and assume that each time series can be reproduced by an epidemic model, two conclusions can be drawn from Figs. 1 and 2: (i) The infected fraction continuously increases (decreases) in the growing (recovering) phases, implying that the initial conditions for each are not randomly chosen but are inherited from the last values of the state variables at the previous , i.e., the succession of status when changing . This is in contrast to most of the previous studies in which the initial conditions for each are always chosen randomly. To reflect the seasonal influence, we may assume that the growing phase corresponds to an adiabatic increase of . For the same reason, we may assume that the recovering phase corresponds to an adiabatic decrease of . (ii) Each burst usually lasts several weeks or months, indicating the time dependence of on a large scale. This is also in contrast to most of the previous studies, where is allowed to change only after the transient process. A key problem is how depends on the time . We notice from our observation that we generally have a few continuously sunny or rainy days in a season, indicating that can be approximately considered as a constant in this short period of time, although changes with on a larger scale. To reflect this feature of real data, we let be a step change, i.e., will be updated to a new value once the system has stayed at the current for a time period , and it can be expressed as

Our model can be any known epidemic model with the above two characteristic features. We take the classical susceptible-infected-susceptible (SIS) model as an example, but a similar analysis can also be performed in the classical susceptible-infected-refractory (SIR) model. In this model, we initially choose a small value of and a few infected seeds. Then, we let the system freely run a period of time where a susceptible individual will become infected with probability if he/she is connected to an infected neighbor, and an infected one will recover to susceptible again with probability  [37]. When there are more than one infected neighbors, a susceptible individual will become infected with probability , where is the number of its infected neighbors. After a time period of , we let have an increase as in Eq. (1), but we keep the individual states unchanged. We repeat this process until reaches its maximum . After that, we simulate the recovering phase by letting but we let the individual states remain. Then, we let the system run as a traditional SIS model. Once for , we let decrease as in Eq. (1), and we keep all the other aspects unchanged. We repeat this process until reaches zero.

NUMERICAL SIMULATIONS

In numerical simulations, we take the random Erdös-Rényi (ER) network with size and average degree as an example, which is generated by the algorithm of Ref. [51]. First, we discuss the case of the SIS model, i.e., we let each node of the network be a model with our two characteristic features. We set and in this paper if there is no specific illustration. Initially, we randomly choose a small fraction of individuals such as and set them as infected while the rest of the population is set as susceptible or healthy.

We increase adiabatically with an increment from to 1, we let the system stay at each for only a finite evolution time , and we compute the infected fraction for each . This is the growing phase. Then, we begin the recovering phase, i.e., we decrease adiabatically with the same and , and we compute for each . Our numerical simulations show that the two phases are generally not overlapping, i.e., with a hysteresis loop. Figure 3(a) shows the results for different , where the curves with squares, circles, up triangles, and down triangles represent the cases of , 2, 5, and 1000, respectively. Comparing the hysteresis loops in Fig. 3(a) with that in Fig. 2(b) we see that they are similar, indicating that the phenomenon of the hysteresis loop has been successfully reproduced by our modified SIS model.

FIG. 3.

(a) and (b) The SIS model. (c) and (d) The SIRS model. (a) and (c) Hysteresis loops for the SIS and SIRS models, respectively, where the curves with squares, circles, up triangles, and down triangles represent the cases of , 2, 5, and 1000, respectively. (b) and (d) and vs for the SIS and SIRS models, respectively. The results are an average of 1000 realizations.

On the other hand, we have noticed that there is a bifurcation point for each hysteresis loop in Fig. 3(a) where the growing and recovering phases will merge into one phase when . The inset of Fig. 3(b) shows the dependence of on , which is a monotonous decreasing curve. We have also noticed that in the case of , the two curves with down triangles in Fig. 3(a) are overlapping, implying no hysteresis loop. To understand this better, we have calculated the hysteresis area surrounded by the growing curve and the recovering curve , i.e., , as shown in Fig. 3(b). It is easy to see that is approximately zero for , indicating that the hysteresis loop exists only when is in the transient process with . By checking the behaviors of , we find that the system has reached its stationary state after . As the stationary state is generally independent of initial conditions, we conclude that the hysteresis loop comes from a memory effect on the adiabatical process with .

Secondly, we discuss the case of the SIR model. Considering that the hysteresis loop exists in different real data of Fig. 1, it is necessary to reproduce it in other epidemic models. For this purpose, we take another classical epidemic model, i.e., the SIR model, as an example. To make the SIR model suitable for the case of recurrent epidemics, we extend it to the susceptible-infected-refractory-susceptible (SIRS) model [37], where an infectious will become a refractory with probability , a refractory will recover to susceptible again with probability , and the other parts will remain the same as in the SIS model. We let and , and we perform the same simulations as in the SIS model. We find that there is also a hysteresis loop in this SIRS model. Figure 3(c) shows the results, where the curves with squares, circles, up triangles, and down triangles represent the cases of , 2, 5, and 1000, respectively. Comparing Fig. 3(c) with Fig. 3(a), we see that they are similar to each other, confirming the universality of the hysteresis loop in both the SIS and SIRS models. We have also calculated the hysteresis area and the bifurcation point for different in the SIRS model, as shown in Fig. 3(d). It is easy to see that there also exists a in Fig. 3(d) where is approximately zero when , confirming again that one of the necessary conditions to observe the hysteresis loop is .

Up to now, the considered growing and recovering phases have been assumed to be separated by the maximum at . An interesting question would be whether there is another way to explain the hysteresis loop in Fig. 2(b). As we mentioned before, may vary from day to day in the real data. Because of the uncertain environment, may not go to before starting the recovering phase. Once this happens, do we still have the hysteresis loop? To answer this question, we let be the turning point to separate the growing and recovering phases. That is, we let increase adiabatically from to and then we let it decrease adiabatically. We first consider the case of the SIS model. Figure 4 shows the numerical simulations, where (a)–(d) represent the cases of , 2, 5, and 10, respectively, and the curves with squares, circles, up triangles, and down triangles in each panel represent the cases of , 0.2, 0.15, and 0.1, respectively. From Fig. 4(a), we see that the case of is similar to Fig. 3(a), implying that this is larger than the of . However, the other three cases of Fig. 4(a) belong to the case of , and they have different behaviors from Fig. 3(a), where their do not reach the maximum at but rather during the recovering phase. That is, the recovering phase is not monotonously decreasing but has a bell shape, indicating a delay effect on the appearance of the maximum . This delay effect can be understood as follows. When is decreased to a value of but still satisfies , the epidemic will keep growing until . The delay effect is common in all four panels of Fig. 4, indicating that it is general for . On the other hand, we notice that the three curves with , 0.15, and 0.1 in Fig. 4(a) keep rising from Fig. 4(b) to Fig. 4(d) and gradually become the normal hysteresis loops. This is because decreases with the increase of [see the insets of Figs. 3(b) and 3(d)] and thus gradually changes from to when changes from Fig. 4(a) to Fig. 4(d).

FIG. 4.

Cases of the SIS model with , where (a)–(d) represent the cases of , 2, 5, and 10, respectively, and the curves with squares, circles, up triangles, and down triangles represent the cases of , 0.2, 0.15, and 0.1, respectively. The results are an average of 1000 realizations.

Then, we discuss the case of the SIRS model. Figure 5 shows the results. Comparing Fig. 5 with Fig. 4, we see that they are similar to each other, confirming that the case of can also show the hysteresis in different epidemic models.

FIG. 5.

Case of the SIRS model, where (a)–(d) represent the cases of , 2, 3, and 10, respectively, and the curves with squares, circles, up triangles, and down triangles represent the cases of , 0.12, 0.16, and 0.2, respectively. The results are an average of 1000 realizations.

All the above results are obtained for fixing . Do they depend on the value of ? To figure out the answer, we change but keep the other parameters unchanged. Figure 6(a) shows the results for the case of the SIS model where the curves with squares, circles, and triangles represent the cases of , 0.001, and 0.0001, respectively. We see that the area of the hysteresis loop decreases with decreasing , indicating that the hysteresis loop will disappear in the limit . We arrived at the same conclusion for the case of the SIRS model; see Fig. 6(b) for the corresponding results.

FIG. 6.

Dependence of the hysteresis loop on the chosen with , , , and , where the curves with squares, circles, and triangles represent the cases of , and 0.0001, respectively. Parts (a) and (b) represent the cases of the SIS and SIRS models, respectively, with in (b).

A BRIEF THEORETICAL ANALYSIS

To gain a deeper insight into the origin of the hysteresis loop in recurrent epidemics, we perform a brief theoretical analysis here by using the SIS model as an example. As the random ER network can be considered homogeneous [37], we take the framework of mean-field theory. For convenience, we divide the evolution time into segments . Let represent the value of at time . For an adiabatic process, we have , with being the increment, and the “” representing the growing and recovering phases, respectively. The evolution of during the segment can be described by the following mean-field equation: By dividing on both sides of Eq. (2), we have Letting be a new variable and integrating Eq. (3) for an evolution time , we obtain its solution as where and is the initial value of for the segment and also the last value of in the segment with . This dependence of on time is in contrast to previous cases in which is a constant for an evolution. On the other hand, the growing and recovering phases have the same range of varying . Do the two phases have the same when they pass the same value of ? The answer is no, and the reason is as follows. In the growing phase, is the minimum of in the segment , while in the recovering phase, is the maximum of in the segment . Therefore, the way to choose in Eqs. (4) and (5) leads to asymmetric initial conditions for the two phases of growing and recovering in the segment with the same , which is how the memory effect occurs, and this is the direct reason why the hysteresis loop is produced. The down triangles and pentagons in Fig. 7 show the solutions of Eq. (4) with , , and , where the down triangles represent the case of , and the pentagons represent the case of . For comparison, we also plot the corresponding numerical results in Fig. 7, where the squares and left triangles represent the cases of and , respectively. It is easy to see that the numerical simulations and theoretical results are qualitatively consistent.

FIG. 7.

Comparison between numerical simulations and theoretical approaches of mean-field theory (MFT), the microscopic Markov-chain approach (MMCA), and effective degree theory (EDT) with , , and , where the squares, circles, up triangles, and down triangles represent the case of for the four approaches, respectively, and the left triangles, right triangles, diamonds, and pentagons represent the case of for the four approaches, respectively. The results for the simulations are an average of 1000 realizations.

For achieve better consistency between the numerical simulations and the theoretical predictions, we have also introduced the mechanism of memory into the microscopic Markov-chain approach (MMCA) [52,53] and the effective degree theory (EDT) [28,53,54] by considering the adiabatical increase of ; see Appendixes A and B for details. The circles, up triangles, right triangles, and diamonds in Fig. 7 show the results of these two methods for and 0.15, respectively. We see that they are more consistent with the numerical simulations than the mean-field theory (MFT) of Eq. (4).

The dynamical equations of both the SIS and SIRS models used in this work are deterministic, indicating that there is no time direction and thus no difference between growing and recovering phases. However, we show analytically that with the introduction of , i.e., the effect of adiabatic continuation, both models become asymmetric on the growing and recovering phases. This asymmetry induces a dependence of dynamics on the integration paths or equivalently a time direction, which may be a new feature in the field of epidemic spreading.

DISCUSSION AND CONCLUSIONS

Both the numerical simulations and the theoretical analysis have confirmed the hysteresis loop in real data. A characteristic feature of this kind of hysteresis loop is that both its growing and recovering phases are smooth, indicating a second-order transition. This feature is different from the previous results in Refs. [44,45], where a first-order transition with a hysteresis loop is predicted in their models. The data we collected in this paper support the features of a second-order transition, but some data may also exist with the features of a first-order transition, which may be found in the near future.

We notice from Fig. 7 that although the effective degree theory is better than the mean-field theory, it is still not completely consistent with the numerical simulations. In fact, this effective degree theory is quite similar to the approach of improved compartmental formalism [55] in which nodes are categorized not only by their state of infectiousness, but also by the state of their neighbors. We also notice that in Ref. [55], the approach works quite well for Poisson networks. This causes us to wonder why the effective degree theory does not perform well at predicting the results of the numerical simulations. To figure out the answer, we apply the effective degree theory to the Poisson networks, and surprisingly it also works well. Thus, we conjecture that the gap between the numerical simulations and the analytical approaches in Fig. 7 comes from the adiabatical feature of evolution. Roughly speaking, all of the theoretical approaches, including both the effective degree theory and the approach in Ref. [55], are based on probability theory, which is good for smooth processes. However, the changing of in an adiabatic way somehow breaks this feature of “smooth” and thus results in the small gap between the theories and simulations in Fig. 7.

In conclusion, we have presented an approach to study the mechanism of recurrent epidemics in real data. We show that there is a hysteresis loop in real data, and it can be reproduced by incorporating two characteristic features into a general epidemic model, i.e., updating before the system is stabilized and adopting the system's state with its previous . Two ways to generate the hysteresis loop are presented: in the first, the growing and recovering phases are separated by a maximum of infected fraction, while in the second, the infected fraction keeps increasing even in the recovering phase. A theoretical analysis shows that the hysteresis loop comes from a memory effect on the adiabatic process with , which enhances our understanding of real epidemics.