Multiple peaks patterns of epidemic spreading in multi-layer networks.

The study of epidemic spreading on populations of networked individuals has seen recently a great deal of significant progresses. A common point in many of past studies is, however, that there is only one peak of infected density in each single epidemic spreading episode. At variance, real data from different cities over the world suggest that, besides a major single peak trait of infected density, a finite probability exists for a pattern made of two (or multiple) peaks. We show that such a latter feature is distinctive of a multilayered network of interactions, and reveal that a two peaks pattern may emerge from different time delays at which the epidemic spreads in between the two layers. Further, we show that the essential ingredient is a weak coupling condition between the layers themselves, while different degree distributions in the two layers are also helpful. Moreover, an edge-based theory is developed which fully explains all numerical results. Our findings may therefore be of significance for protecting secondary disasters of epidemics, which are definitely undesired in real life.

Introduction

Epidemic spreading in networked populations has been studied intensely in the last decade, and a lot of great progresses has been achieved [1], [2], [3], [4] which significantly increased our understanding. This is actually useful for public health authorities to assess situations quickly, to take and enforce informed decisions, and to optimize vaccination and drug delivery policies. Initially, the main attention focused on static networks, where each node represents an immobile agent and the contagion occurs only between neighboring nodes: it was revealed that scale-free networks display a vanishingly small epidemic threshold in the thermodynamic limit [5], [6]. After that, the focus shifted to reaction-diffusion models [7], [8], [9], [10], flow-driven epidemics [11], [12], [13], [14], [15], objective spreading [16], [17] and adaptive behaviors [18], [19], [20], [21], [22]. Finally, in a third stage, multilayered [23], [24], [25], [26], [27], [28], [29], [30], [31], [32] and temporal [33], [34] networks were assumed to play a critical role on such processes.

A common feature in past studies is the use of phenomenological models which typically produce a single peak of infected density in each individual epidemic spreading. An interesting question is therefore whether or not all real evolutionary processes are conveniently represented by such a framework. A scrupulous analysis of a large number of real data from different cities over the world surprisingly shows that, besides a major pattern made of a single peak, there is a finite and non negligible probability for a new pattern of epidemic outbreak featuring two (or multiple) peaks. Notice that a two peaks pattern implies two outbreaks in a single spreading period, i.e. a secondary occurrence of the same epidemics, which may in turn produce severe calamities and disasters within unprepared populations. Understanding the underlying mechanism at the basis of this new pattern (with the help of a novel model extracted from real data) is therefore quintessential to properly cope with such life-threatening hazards. Previous studies showed that the phenomenon of multiple peaks in epidemics may occur in weak coupled different communities, where a full homogeneous mixing cannot be assumed [35], [36], [37], [38]. Considering that human behaviors occur in social networks, a key problem for further studies is how the network structures influence the formation of two-peaks patterns and the transition between the case of single peak and that of two peaks.

We here study this problem by presenting a network based model to reproduce the phenomenon of two peaks from epidemic data. As social activities and interactions occur in network structures, we here consider the epidemics in different geographic regions (or cities) as that occurring in multilayered graphs. Namely, we will take two coupled neighboring regions as an example, and construct a two-layered network model which fully reproduces the observed patterns of epidemics. In particular, we demonstrate that the pattern of two peaks originates from large time delays of epidemic outbreaks between the two layers, which depends in turn on both the difference in the degree distributions of the two layers and a weak coupling condition between them. To better understand the findings, an edge-based theory is developed which perfectly agree with the numerical simulations.

Patterns of epidemic outbreak with two or multiple peaks in real data

Monitoring the potential outbreaks of an epidemic spreading is of extreme importance for protection of our society. Based on the detected trend of spreading of infections such as SARS (Severe Acute Respiratory Syndrome), H1N1 (Swine Influenza), H5H1 (Avian Influenza), and Ebola, one can indeed attempt to enforce suitable measures able to reduce the epidemic at its maximum extent. For this purpose, many countries have established their sentinel surveillance systems to collect epidemic data. For instance, Hong Kong Department of Health has organized a surveillance system, with the aim of collecting empirical data of infectious diseases, and of analyzing and predicting the trend of the infection. In such a system, there are about 64 General Out-Patient Clinics (GOPC) and 50 General Practitioners (GP), which form two distinct sentinel surveillance networks of the city [39], [40]. In these two networks, one obtains for instance the weekly consultation rates of influenza-like illness (per 1,000 consultations), which reflect the overall influenza-like illness activity in Hong Kong.

Fig. 1 (a) shows the collected data from 1998/1/3 to 2014/8/2 in GP, while the corresponding data of GOPC is not shown here. From Fig. 1(a) one sees that there are many events of epidemic spreading, and the intervals between two consecutive events are not regular, indicating non-periodic outbreaks of recurrent epidemics [40]. On the other hand, one notices from Fig. 1(a) that most of the outbreaks correspond to a single peak of infected density, which is the pattern well described by the classical susceptible-infected-refractory (SIR) models. However, one notices also that there is a finite probability for a novel pattern of epidemic outbreak which features, instead, two or multiple peaks. Fig. 1(b) shows one of such patterns (with two peaks) occurring at around 2005/6, which indicates that an infectious disease raised two times during that epidemic period in Hong Kong. Such unexpected phenomenon also exists in the data from GOPC (not shown here).

Fig. 1

Time series of recurrent epidemics in different cities over the world. (a) The weekly consultation rates of influenza-like illness (per 1000 consultations) from 1998/1/3 to 2014/8/2 in Hong Kong for the General Practitioners (GP) sentinel system. (b) Zoom of one of the patterns with two peaks, occurring at around 2005/6 in (a). (c) The time series of reported weekly measles infective cases I in Boston. (d) Zoom of one of the patterns with two peaks, occurring at around 1915/4 in (c). (e)–(h): Time series of infectious disease in two coupled cities. (e) and (f): The yellow line represents the total number of weekly measles infective cases in coupled Bristol and Newcastle, while the blue and green lines represent the number in Bristol and Newcastle, respectively. (f) is one of the patterns, occurring at around 1957/3 in (e). (g) and (h): The yellow line represents the total number of weekly measles infective cases in coupled Bristol and Sheffield, while the blue and green lines represent the number in Bristol and Sheffield, respectively. (h) is one of the patterns, occurring at around 1955/5 in (g). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

These multiple peaks patterns are actually occurring generically, and are not limited to a specific geographical region. Epidemic data from other sources and cities display, indeed, ubiquitously patterns similar to that reported in Fig. 1(b). For instance, Fig. 1(c) shows the data of the weekly measles infective cases (WMICs) from 1908 to 1937 in Boston [41], [42], and a two peaks pattern is shown in Fig. 1(d) at around 1915/4. In addition, two or multiple peaks of infected density may characterize the outbreaks in the total number of WMICs of two neighboring cities, also when uni-modal patterns are actually observed in each individual city. For example, Fig. 1(e) and (f) report the data of WMICs in Bristol and Newcastle [43]: the yellow line denotes the total number of WMICs, whereas the blue and green lines represent the data in Bristol and Newcastle, respectively. In Fig. 1(f) one can well appreciate that the pattern of two peaks occurs only in the total number of WMICs. Similarly, Fig. 1(g) and (h) show the case of Bristol and Sheffield [43], and once again a typical two peaks pattern [Fig. 1(h)] occurs.

The two-layered network model

To capture the underlying mechanism, we introduce a model of a two-layered network, where the two layers represent actually two interconnected regions or cities. Fig. 2 is a sketch of the model: and are the two layers, which are coupled through the inter-network . For the sake of simplicity, we let the two networks and have the same size . Furthermore, ⟨k⟩, ⟨k⟩, and ⟨k⟩ represent the average degrees of networks and respectively. In details, we first generate two separated networks and with the same size N and different degree distributions P(k) and P(k), respectively. Then, we add links between and . That is, we randomly choose two nodes from and and then connect them if they are not connected yet. The process is repeated until all the needed specifications are attained. In the above way, one obtains a uncorrelated two-layered network.

Fig. 2

Sketch of the two-layered network model, which reproduces the pattern of two peaks. “Black”, “green” and “red” lines represent the links of the networks and the inter-network respectively. β and β denote the infectious rates of networks and . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Under such a framework, we adopt the Susceptible-Infected-Refractory (SIR) model to study the epidemic spreading dynamics. Particularly, each node is a unit of the SIR model, where S, I and R represent the susceptible, infected and refractory phases of individuals, respectively. At each time step, a susceptible node will be infected by an infected neighbor with rate β, and an infected node will become refractory with probability μ. The infectious process will be considered terminated when no more infected nodes exist. While μ is taken to be the same for all networks, we let β and β be the infectious rates of networks and respectively. Thus, a susceptible node will be infected with a probability where k is the infected neighbors in the same network and is the infected neighbors in the coupled network.

In our numerical simulations, we choose to be a scale-free (SF) network with degree distribution [44], and a random regular (RR) network with a constant degree k [45]. The network size is the average degree and initially 0.1% of individuals of network are chosen to be infected.

Results

The two peaks outbreak pattern in simulations

A pattern of two peaks appears in the numerical simulations of our model. The yellow circles in Fig. 3(a) report the evolution of the infected density ρ in the whole network with the parameters and . One can easily differentiate two peaks of ρ, indicating that the empirical observations in Fig. 1 can be fully reproduced.

Fig. 3

The two peaks pattern occurring in different conditions, with and where the symbols represent the simulated results and the lines denote the corresponding theoretical predictions (calculated via the edge-based compartmental theory described in the Method section). (a) ρ(t) vs. t, where ρ, and represent the infected densities in the entire network, and the networks and respectively. Other parameters are and . (b) The influence of ⟨k⟩ on ρ(t) with and where the “triangles”, “squares” and “circles” represent the cases of and 0.5, respectively. (c) The influence of β on ρ(t) with and where the “triangles”, “squares” and “circles” represent the cases of and 0.001, respectively. (d) The influence of γ on ρ(t) with and where the “triangles”, “squares” and “circles” represent the cases of and 1.5, respectively.

To gather a deeper understanding, we also measure the evolutions of and in both layers and and report them as green triangles and blue squares in Fig. 3 (a), respectively. It is easy to see that the times at which the maximum infected density is obtained in layers and are different, implying that the pattern of two peaks is likely triggered by the time difference of epidemic outbreaks in the two layers.

The factors determining the occurrence of two peaks pattern

The next step is to focus on the factors that determines the occurrence of such two peaks pattern. For this purpose, we first concentrate on the role of the average degree ⟨k⟩ of the inter-network . Fig. 3(b) shows the evolution of the infected density ρ with different ⟨k⟩ (with triangles, squares and circles denoting the cases of and 0.5, respectively). One notices that the pattern of ρ is uni-modal when ⟨k⟩ is large, but bimodal when ⟨k⟩ is sufficiently small, indicating that ⟨k⟩ is a key factor for the appearance of a bimodal pattern: a smaller ⟨k⟩ favours the appearance of the two peaks pattern, indicating that the two main networks and should be only weakly coupled among them.

Further, we study the influence of the infectious rate β on the pattern. Fig. 3(c) reports the results obtained for different β (with triangles, squares and circles representing the cases of and 0.001, respectively). Once again one may notice that the condition of a weak coupling is essential for a bimodal pattern: ρ is indeed uni-modal when β is large, while the two peaks appear when β is small. Finally, we study the influence of the exponent γ of the SF network. Fig. 3(d) shows that the bimodal feature is reduced with the increase of γ. As a larger γ means a smaller difference between the structures of the SF and RR networks, one can infer that the heterogeneity between the two layers is very helpful for the appearance of two peaks patterns.

All the numerical results are fully confirmed by an edge-based compartmental theory, see our theoretical Eqs. (18) and (19) in the Theoretical analysis section. The solid curves in Fig. 3(a)–(d) show the corresponding theoretical results.

These numerics point that a weak coupling is the necessary condition for the emergence of the new pattern, while a difference in heterogeneity between the two layers is a very helpful factor. An interesting question is whether the observed behavior corresponds to a critical phenomenon. To figure out the answer, we let τ be the time interval between the two peaks of ρ. In particular, the two peaks will merge into a single one when . Similarly, we let be the time delay between the two peaks in layers and where and are the times of occurrence of the peak in and respectively. The trivial situation would be that for which but our numerical simulations show that this condition is attained only when τ is large, whereas one has τ < δ when τ is sufficiently small. And, in particular, one has δ > 0 when .

Fig. 4 (a) and (b) show the dependence of τ and δ on β for fixed and different ⟨k⟩, respectively. From Fig. 4(a), one sees that when ⟨k⟩ is small, τ will decrease monotonically and be non-vanishing with the increasing of β, indicating that the event of two peaks always exists in the pattern. However, when ⟨k⟩ is increased, τ will decrease rapidly to zero, implying that there is a critical for different ⟨k⟩. When the epidemic spreading event occurs through a pattern of two peaks in the infected density, while for it occurs via the traditional pattern with a single peak. On its turn, Fig. 4(b) shows that δ decreases monotonically with the increase of β for all the three cases of ⟨k⟩, and it never vanishes. This is because the spreading speed is different in homogeneous and heterogeneous networks. Generally speaking, epidemic spreading is faster in heterogeneous network than in homogeneous network [46].

Fig. 4

(a) and (b) report the dependence of τ and δ on β (see main text for definitions) with different ⟨k⟩ in SF-RR networks. (c) and (b) report the dependence of τ and δ on β with different γ in SF-RR networks. Symbols represent the simulated results and the lines are the corresponding theoretical results (calculated via the edge-based compartmental theory described in the Method section). Other parameters as in the caption of Fig. 3. All the results are averaged over 100 independent realizations.

We then turn to investigate the influence of the heterogeneity in degree distribution on the occurrence of the two peaks pattern. Fig. 4(c) and (d) show the dependence of τ and δ on β for different γ and fixed . While the network heterogeneity decreases with the increasing of γ, it is easy to see that when γ is large, τ and δ decrease more prominently with β. Specifically, the difference of spreading speeds between the two layers is not distinctive for large γ, resulting in the disappearance of the two peaks. Therefore, increasing the coupling strength (i.e. ⟨k⟩ and β) or decreasing the heterogeneity of network topology between the networks and will decrease the time delay of epidemic outbreak and then suppress the pattern of two peaks.

We have also confirmed all these numerical results in Fig. 4(a)–(d) by the theoretical Eqs. (18) and (19) in the Theoretical analysis section. For each set of parameters in Fig. 4(a)–(d), we first produce the infected densities and from Eqs. (18) and (19), as done in Fig. 3(a)–(d), and then measure the corresponding τ and δ. The solid curves in Fig. 4(a)–(d) show the theoretical τ and δ. One can easily see that the theoretical results are fully consistent with the numerical results.

The effect of the coupling strength on the final epidemic size

In the above text, we have discussed that the value of the coupling strength is crucial for determining the kind of pattern of the epidemic spreading. Then, a natural question is how the coupling strength influences the final size of epidemic spreading. Fig. 5 shows the dependence of the final value of ρ on β and β for a scale-free (SF)-random regular (RR) network configuration, with different ⟨k⟩ and β. All the results are averaged over 100 independent realizations with and in (a), and in (c), and and in (e), respectively. Comparing Fig. 5(a) with (c), one sees that the blue area in (c) is much smaller than that in (a), indicating that increasing β has the effect of enhancing the final epidemic size. Comparing Fig. 5(a) with (e), one sees that the blue area in (e) is slightly smaller than that in (a), indicating that increasing ⟨k⟩ can also slightly enhance ρ.

Fig. 5

Dependence of the final epidemic size ρ on β and β in SF-RR networks with and in (a)-(b), and in (c)-(d), and and in (e)-(f). Left and right panels correspond to numerical simulations and theoretical results, respectively. All the numerical results are averaged over 100 independent realizations. Other parameters are set as and .

On the other hand, one can also obtain ρ by the theoretical formulae in the theoretical analysis Section. The right panels in Fig. 5 (i.e. Fig. 5(b), (d) and (f)) show the results corresponding to the left panels, respectively, and indicate that all theoretical predictions are highly consistent.

Theoretical analysis

Edge-based compartmental theory for a single network

Let us first illustrate the edge-based compartmental theory for a single network, by following the methods and tools introduced in Refs. [47], [48], [49], [50], [51], [52], [53], [54], [55], [56].

For an uncorrelated, large and sparse network, the SIR model can be described in terms of the quantities S(t), I(t), and R(t), which represent the densities of the susceptible, infected, and recovered nodes at time t, respectively. Let θ(t) be the probability that a neighbor v of u has not transmitted the disease to u along the edge connecting them up to time t. Then, the node u with degree k is susceptible at time t as . Averaging over all k, the density of susceptible nodes at time t is given by where P(k) is the degree distribution of the network. In order to solve for S(t), one needs to know θ(t). Since a neighbor v of node u may be susceptible, infected, or recovered, θ(t) can be expressed as where Φ(t), Φ(t), Φ(t) is the probability that the neighbor v is in the susceptible, infected, recovery state, respectively, and has not transmitted the disease to node u through their connection. Once these three parameters can be derived, we will get the density of susceptible nodes at time t by substituting them into Eq. (2) and then into Eq. (1). To this purpose, in the following, we will focus on how to solve them.

To find Φ(t), we now consider a randomly chosen node u, and assume this node is in the cavity state, which means that it cannot transmit any disease to its neighbors v but can be infected by its neighbors. In this case, the neighbor v can only get the disease from its other neighbors except the node u. Thus, node v with degree k′ is susceptible with probability at time t. For uncorrelated networks, the probability that one edge from node u connects with a node v with degree k′ is k′P(k′)/⟨k⟩. Summing over all possible k′, one obtains

According to the SIR spreading process, the growth of Φ(t) includes two consecutive events: first, an infected neighbor has not transmitted the infection to node u via with probability ; second, the infected neighbor has been recovered with probability μ. Combining these two events, the Φ(t) to Φ(t) flux is . Thus, one gets

Once the infected neighbor v transmits the disease to u successfully (with probability β), the Φ(t) to flux will be βΦ(t), which means That is Combining Eqs. (4) and (5), and considering (as initial conditions) and one obtains Substituting Eqs. (3) and (6) into Eq. (2), one gets an expression for Φ(t) in terms of θ(t), and then one can rewrite Eq. (5) as

After obtaining the solution of Eq. (7), we can substitute it into Eq. (1) to obtain S(t). Then, I(t) and R(t) can be obtained as follows In fact, Eq. (7) does not depend on Eq. (8), so the system is governed by the single ordinary differential equation (7). Although the resulting equation are simpler than those found by other methods, it can be proven to exactly predict the disease dynamics in the large-population limit for different network topologies [53], [57].

The theory for two-layered networks

When one assumes that the population is made up of two interacting networks, then P(k 1, k 2) denote the probability that a node of network j has k 1 degree in network 1 and k 2 in network 2. For the sake of simplicity, one can name the two networks and as 1 and 2. Let β be the rate of transmission across an edge from network l to network j, and let us define μ to be the recovery rate of a node in any network.

θ can be defined to be the probability that an edge to a test node u in network j () is coming from network l (), and has not transmitted the infection.

Now, θ 1, 2 can be solved as in the case of a single network. Since a neighbor v in network 2 of node u in network 1 may be susceptible, infected, or recovered, θ 1, 2 can be expressed as where is the probability that the neighbor v is in the susceptible, infected, recovery state, and has not transmitted the disease to node u through their connection.

Similarly, to find the neighbor v in network 1 can only get the disease from its other neighbors except the node u in network 2. Thus, the node v with degree k 1 in network 1 and degree k 2 in network 2 is susceptible with probability at time t. For uncorrelated networks, the probability that one edge from node u connects with a node v with degree (k 1, k 2) is . Thus, one has

It is easily to know that the growth of includes two consecutive events: first, an infected neighbor has not transmitted the infection to node u via with probability ; second, the infected neighbor has been recovered with probability μ. Combining these two events, the to flux is . Thus, one gets

Once the infected neighbor v in network 1 transmits the disease to node u in network 2 successfully (with probability β 1, 2), the to flux will be which means

Combining Eqs. (11) and (12), and considering (as initial conditions) and one obtains So, one gets

Similarly, one can write down θ 1, 1, θ 2, 1 and θ 2, 2 as follows With Eqs. (14)–(17) on hand, the densities associated with each distinct state can be obtained by

Eqs. (18) and (19) are the main theoretical results from which the theoretical results in Fig. 3, Fig. 4, Fig. 5 are calculated. Now, a critical question comes out: can our theory find out the epidemic threshold? In fact, we find that the threshold for the whole network to show epidemic outbreak can be theoretically figured out by the Jacobian matrix of Eqs. (14)–(17). Particularly, the epidemic thresholds and can be obtained by the linear stability analysis of the trivial fixed point solution of Eqs. (14)–(17). Without loss of generality, one can set the recovery rate . The Jacobian matrix of Eqs. (14)–(17) at (1, 1, 1, 1) can be expressed as

For uncorrelated networks, the matrix can be simplified as and can be obtained by seeking for the condition where the largest eigenvalue Λ of changing from 0 to positive. Instead of obtaining the expression of Λ of we here numerically find the transition point of Λ from 0 to positive. The white lines in Fig. 6 show the results on SF-RR networks, i.e. the thresholds and where (a) and (b) represent the case of and (c) and (d) the case of ; and (a) and (c) represent the case of and (b) and (d) the case of . For confirmation of such a prediction, one can use the same set of parameters, and calculate numerically ρ of the whole SF-RR networks, which is shown in Fig. 6(a)–(d). The phase diagram of ρ in each panel of Fig. 6(a)–(d) can be divided into two distinct regions (whose boundaries are depicted in white lines): the healthy and endemic phases, respectively. The theoretical predictions are in good agreement with the numerical results.

Fig. 6

Phase diagram of the final epidemic size. Dependence of ρ on β and β on SF-RR networks with in (a)-(b) and in (c)-(d). Left and right panels correspond to and respectively. The Color-coded values ρ are obtained from numerical simulations and the theoretical predictions (white line) are the solutions of the largest eigenvalue of Jacobian matrix in matrix (20). All numerical results are averaged over 1,000 independent realizations. Other parameters are and .

When networks and can be treated as isolated graphs (i.e. when their coupling is particularly weak), the explicitly expression of the epidemic thresholds is obtained by setting i.e. and . When the coupling between the layers is increased, and decrease gradually. This result is consistent with the evidences reported in Ref. [25], [58].

Fig. 7 shows the phase diagram of epidemic for different ⟨k⟩ and β, where (from top to bottom) the curves with different colors represent the epidemic threshold lines of and for the cases of respectively. Panels (a) and (b) correspond to and 2.0, respectively. In detail, each of the curves represents the boundary between the healthy and endemic phases, and the solutions are obtained by solving the largest eigenvalue of the Jacobian matrix . Other parameters are set as and . Once the coupling strength is strong enough, the value of Λ is mainly determined from the contribution of the interconnection network which is similar to a pure bipartite network. Looking at Fig. 7, an interesting conclusion is that an endemic state can occur even if there are no interactions within each individual network. These results imply that, in order to control an epidemic spreading, one needs to pay more attention to the effects of the coupling strength between different regions.

Fig. 7

Phase diagrams of the healthy and the endemic phase for (a) and (b) . The curves indicate the boundaries of healthy and the endemic phases, where the solution are obtained by solving in matrix (20). Different color curves from top to down correspond to the case of . Other parameters are and .

Discussion

Let us remark that the weak coupling condition predicted by us for the occurrence of the novel epidemic pattern is actually consistent with the cases of the data of Fig. 1. As it is well known, indeed, Hong Kong in Fig. 1(a) consists of islands, and the movement of individuals between different islands is not as convenient as that within each single island, and thus the coupling between neighboring islands can be considered to be weak. At the same time, the population distribution in Hong Kong central island is significantly different from that characterizing the surrounding islands, confirming the presence of the second ingredient predicted by our theory: i.e. a difference in the heterogeneity of the layers’ structures.

In Boston, a river separates the city into two parts, which (to all extent) can be considered as equivalent to two islands. The same reasoning applies to the neighboring cities of Bristol and Newcastle and the neighboring cities of Bristol and Sheffield. As Bristol and Newcastle, Newcastle and Sheffield are separated regions in the United Kingdom, and they can therefore be considered as a pair of weakly coupled networks.

Our predictions were obtained on coupled SF-RR networks, and it is legitimate to seek for generality of the two peaks pattern phenomenon, by means of investigating coupled networks with other topological structures. For this purpose, we have also studied the case of SF-SF and RR-RR networks, respectively. Very interestingly, one finds that the pattern of two peaks can be still observed by adjusting the coupling strength between the coupled layers (see Fig. 8 ). Although the two coupled SF networks (or the two coupled RR networks) can be considered as the case with same heterogeneity, they are produced independently and thus their network topologies are different, indicating that there exists still the effect of heterogeneity between the two coupled networks. On the other hand, extension to three-layered model was also considered, and it was found that there is a small probability to produce a pattern of three peaks (see Fig. 9 ). Therefore, while in principle one can expect a multi-peaks pattern to occur in a multilayered network, the majority of unusual cases (i.e. cases in which the epidemic event is not happening with a single maximum of infected density) will be characterized by just two peaks, in full consistency with the data of Fig. 1.

Fig. 8

The patterns of two peaks for different layer structures. (a), (b) and (c) represent the time evolutions of the infected density ρ in RR-SF networks, SF-SF networks, and RR-RR networks, respectively. (d), (e) and (f) show the time evolutions of and corresponding to (a), (b) and (c). The parameters in the left, middle, and right panels are and ; and ; and and respectively. Other parameters are and .

Fig. 9

Three peaks patterns in three-layered networks. Panels (a), (b) and (c) report the patterns with one peak, two peaks and three peaks of ρ, for a three-layered SF-SF-RR network. (d), (e) and (f) show the time evolutions of and corresponding to (a), (b) and (c). The parameters in the left, middle, and right panels are ; ; and respectively. Other parameters are and .

In summary, epidemic spreading has been well studied in the past decades but mainly focused on outbreaks corresponding to patterns with a single peak of infected density. We here reported (from real data) the evidence that also a pattern of two peaks in a single epidemic period is possible. We pointed out that such a pattern is a genuine product of a multi-layered interaction structure, and we have introduced a proper model able to fully capture the mechanisms for its occurrence. Our model, together with reproducing the classical pattern of a single peak, can generate the pattern with two peaks when coupling between the layers is weak, while the difference in heterogeneity of the layers’ structures is also helpful.