Modeling correlated bursts by the bursty-get-burstier mechanism
MModeling correlated bursts by the bursty-get-burstier mechanism
Hang-Hyun Jo ∗ Asia Pacific Center for Theoretical Physics, Pohang 37673, Republic of KoreaDepartment of Physics, Pohang University of Science and Technology, Pohang 37673, Republic of Korea andDepartment of Computer Science, Aalto University, Espoo FI-00076, Finland (Dated: March 2, 2018)Temporal correlations of time series or event sequences in natural and social phenomena havebeen characterized by power-law decaying autocorrelation functions with decaying exponent γ . Suchtemporal correlations can be understood in terms of power-law distributed interevent times withexponent α , and/or correlations between interevent times. The latter, often called correlated bursts,has recently been studied by measuring power-law distributed bursty trains with exponent β . Ascaling relation between α and γ has been established for the uncorrelated interevent times, whilelittle is known about the effects of correlated interevent times on temporal correlations. In order tostudy these effects, we devise the bursty-get-burstier model for correlated bursts, by which one cantune the degree of correlations between interevent times, while keeping the same interevent timedistribution. We numerically find that sufficiently strong correlations between interevent times couldviolate the scaling relation between α and γ for the uncorrelated case. A non-trivial dependence of γ on β is also found for some range of α . The implication of our results is discussed in terms of thehierarchical organization of bursty trains at various timescales. I. INTRODUCTION
A variety of dynamical processes observed in nat-ural and social phenomena are known to show non-Poissonian or inhomogeneous temporal patterns. Exam-ples include solar flares [1], earthquakes [2, 3], neuronalfirings [4], and human communication and locomotor ac-tivities [5, 6]. Temporal correlations in such time seriesor event sequences have often been described in termsof 1 /f noise [7–9] or power-law decaying autocorrelationfunctions [10, 11]. The autocorrelation function for anevent sequence x ( t ) is defined with delay time t d as fol-lows: A ( t d ) = (cid:104) x ( t ) x ( t + t d ) (cid:105) t − (cid:104) x ( t ) (cid:105) t (cid:104) x ( t ) (cid:105) t − (cid:104) x ( t ) (cid:105) t , (1)where (cid:104)·(cid:105) t means a time average. The event sequence x ( t )can be considered to have the value of 1 at the momentof event occurred, 0 otherwise. For the event sequenceswith long-term memory effects, one may find a power-lawdecaying behavior with decaying exponent γ : A ( t d ) ∼ t − γd . (2)The autocorrelation function captures the entire tempo-ral correlations present in the event sequence, which canbe understood in terms of interevent times and correla-tions between interevent times. Here the interevent time,denoted by τ , is defined as a time interval between twoconsecutive events.The heterogeneous properties of interevent times haveoften been characterized by the power-law intereventtime distribution P ( τ ) with power-law exponent α : P ( τ ) ∼ τ − α , (3) ∗ [email protected] which may readily imply clustered short interevent timeseven without correlations between interevent times. Thisphenomenon has been described in terms of bursts, i.e.,rapidly occurring events within short time periods al-ternating with long inactive periods [5]. It is well-known that bursty interactions between individuals havea strong influence on the dynamical processes takingplace in a network of individuals, such as spreading ordiffusion [12–17]. In case when interevent times are fullyuncorrelated, i.e., for renewal processes, the power spec-tral density was analytically calculated from power-lawinterevent time distributions [18]. Using this result, onecan straightforwardly derive the scaling relation between α and γ : α + γ = 2 for 1 < α ≤ α − γ = 2 for 2 < α ≤
3. (4)This relation was also derived in the context of priorityqueueing models [19]. In addition, Abe and Suzuki [20]derived α + p = 2 for 1 < α < < p <
1, wherethe scaling exponent p characterizes Omori law in earth-quakes, with a possible interpretation of p = γ .Then a natural question arises: If the interevent timesare correlated with each other, would the above scalingrelation still hold? In other words, in order to violate theabove scaling relation, how strong correlations betweeninterevent times should be introduced? In order to studythis question, the correlations between interevent timescan be measured using the notion of bursty trains [10].A bursty train is defined as a set of events such that in-terevent times between any two consecutive events in thebursty train are less than or equal to a given time window∆ t , while those between events in different bursty trainsare larger than ∆ t . The number of events in the burstytrain is called burst size, and it is denoted by b . Thedistribution of b follows an exponential function if theinterevent times are fully uncorrelated with each other. a r X i v : . [ phy s i c s . d a t a - a n ] S e p However, b has been empirically found to be power-lawdistributed, i.e., P ∆ t ( b ) ∼ b − β (5)for a wide range of ∆ t , e.g., in earthquakes, neuronal ac-tivities, and human communication patterns [10, 21, 22].This indicates the presence of correlations between in-terevent times, hence this phenomenon is called corre-lated bursts . We also note that the exponential distri-butions of P ∆ t ( b ) have recently been reported for mobilephone calls of individual users in another work [23]. Insum, we expect that the temporal correlations observedin the autocorrelation function A ( t d ) can be fully under-stood in terms of the statistical properties of intereventtimes, e.g., P ( τ ), together with those of the correlationsbetween interevent times, e.g., P ∆ t ( b ). Simply put, onecan study the dependence of γ on α and β , which is theaim of this paper.The dependence of γ on α and β has been investigatedby dynamically generating event sequences showing tem-poral correlations, described by the power-law distri-butions of interevent times and burst sizes in Eqs. (3)and (5). These generative approaches were based on ei-ther two-state Markov chain [10] or self-exciting pointprocesses [24]. Although such approaches have been suc-cessful for reproducing empirically-observed correlatedbursts to some extent, our understanding of the under-lying mechanisms for correlated bursts is far from com-plete. In addition to the generative modeling approaches,here we take an alternative approach by devising thebursty-get-burstier model, where power-law distributionsof interevent times and burst sizes are inputs rather thanoutputs of the model. In our model, one can explicitlytune the degree of correlations between interevent timesto test whether the scaling relation in Eq. (4) will be vi-olated due to the correlations between interevent times.Our paper is organized as follows: In Sect. II, thebursty-get-burstier model is devised to simulate eventsequences with tunable correlations between intereventtimes, and to study the effects of such correlations onthe scaling relation established for the uncorrelated in-terevent times. We conclude our work in Sect. III. II. METHODS AND RESULTSA. Uncorrelated interevent times
As a baseline, we investigate the case with uncorrelatedinterevent times. We first relate power-law exponentscharacterizing bursty time signals as mentioned in Sect. I.The power spectral density P ( f ) of a time signal x ( t ),where f denotes the frequency, is known to be the Fouriertransform of the autocorrelation function A ( t d ). Thus if A ( t d ) ∼ t − γd for 0 < γ <
1, then one finds the scaling of P ( f ) ∼ f − η with η = 1 − γ. (6) When the interevent times are i.i.d. random variableswith P ( τ ) ∼ τ − α , implying no correlations between in-terevent times, the power-law exponent η is obtained asa function of α as follows [18, 25]: η = (cid:26) α − < α ≤ − α for 2 < α ≤
3. (7)Combining Eqs. (6) and (7), we obtain Eq. (4): γ de-creases with α for 1 < α ≤ α for 2 < α ≤
3. These power-law exponents can also berelated via Hurst exponent H , i.e., γ = 2 − H [26, 27]or η = 2 H − x ( t ) with n + 1 events using n uncor-related interevent times, { τ , · · · , τ n } , that are indepen-dently drawn from P ( τ ). Let us assume that the zerothevent occurs at t = 0. Then the timing of the i th eventfor i = 1 , · · · , n is given by t i = (cid:80) ii (cid:48) =1 τ i (cid:48) , leading to theevent sequence x ( t ) in the range of 0 ≤ t ≤ t n as x ( t ) = (cid:26) t ∈ { t i } ,0 otherwise. (8)We adopt the power-law interevent time distributionwith α > P ( τ ) = (cid:26) ( α − τ α − τ − α for τ ≥ τ ,0 otherwise, (9)where τ is the lower bound of interevent times. Foreach value of α , 50 event sequences are generated with n = 5 · and τ = 10 − [31]. As shown in Fig. 1, wefind that the numerical value of γ as a function of α iscomparable to the scaling relation in Eq. (4) for α (cid:46) . α (cid:38) .
4. The discrepancy between the theoreticaland numerical values for 1 . (cid:46) α (cid:46) . α = 2 aswell as due to finite-size effects. The finite-size effects canbe studied by measuring the n -dependence of γ as shownin Fig. 1(c), where γ ( n ) = γ ( ∞ ) + cn − with constant c is used to estimate γ ( ∞ ) = 0 . α = 1 . γ ( ∞ ) = 0 . α = 2 .
2, respectively. Both limitingvalues of γ ( ∞ ) are yet different from 0 . α =2 [see Fig. 1(d)], we use the form of γ ( n ) = γ ( ∞ ) + cn − . to estimate γ ( ∞ ) = 0 . n = 5 · for the rest of our paperbecause this number is already large enough to studythe temporal properties of most empirical datasets. In γ α numerictheory 0.28 0.3 0.32 0.34 0.36 0 1x10 -5 -5 (c) γ ( n ) α =1.82.2 0.21 0.24 0.27 0.3 0.33 0 0.01 0.02(d) γ ( n ) α =210 -4 -3 -2 -1 -6 -5 -4 -3 -2 -1 (a) A ( t d ) t d α =1.51.72.02.32.5 FIG. 1. (Color online) Numerical test of the scaling relation inEq. (4) for the event sequences constructed using uncorrelatedinterevent times, that are drawn from P ( τ ) in Eq. (9). (a)Autocorrelation functions A ( t d ) for various values of α . Foreach value of α , we have generated 50 event sequences with n = 5 · and τ = 10 − . (b) The estimated value of γ asa function of α using A ( t d ) ∼ t − γd (red circles), compared tothe theory in Eq. (4) (solid line). (c) Finite-size effects on thevalues of γ ( n ) for α = 1 . .
2, fitted with the form of1 /n . (d) The values of γ ( n ) for α = 2 are better described bythe form of 1 /n . than by that of 1 /n . addition, for all cases, we obtain the exponential burstsize distribution P ∆ t ( b ) ∝ e − b/b c (∆ t ) with exponentialcutoff b c (∆ t ) for the wide range of ∆ t (not shown), asexpected.The decreasing and then increasing behavior of γ asa function of α can be roughly understood by the factthat the autocorrelation function essentially measures thechance of finding two events occurred in t and t + t d , nomatter how many events occur between those two events.This chance can be written as [18] (cid:104) x ( t ) x ( t + t d ) (cid:105) t ∝ ∞ (cid:88) k =1 P (cid:63)k ( t d ) ≡ G ( t d ) , (10)where P (cid:63)k ( t d ) denotes the probability of finding k con-secutive interevent times whose sum is exactly equal to t d . The term with k = 1 in Eq. (10) simply correspondsto the value of the interevent time distribution at τ = t d .This value is proportional to ( α − τ t d ) α that is increas-ing and then decreasing for t d > τ as α varies from 1to 3. Such non-monotonic behavior is expected to otherterms with k > G ( t d ) can have the maximumvalue for an intermediate range of α . With an additionalassumption that the larger values of G ( t d ) for a widerange of t d imply the smaller γ , one can get some hints on why γ has the lowest value around at α = 2. B. Correlated interevent times
In order to implement the correlations between in-terevent times, we devise the bursty-get-burstier modelof correlated interevent times, where the power-law dis-tributions of interevent times and burst sizes are inputsrather than outputs of the model. Along with the power-law form of P ( τ ) in Eq. (9), we adopt power-law burstsize distributions as P ∆ t ( b ) ∼ b − β for a wide range of∆ t , meaning that an event sequence will be constructedfrom the prescribed power-law distributions of intereventtimes and burst sizes. As in the uncorrelated case, weprepare the set of n interevent times, T ≡ { τ , · · · , τ n } ,that are independently drawn from P ( τ ). Correlationsbetween interevent times are implemented by shufflingor permuting those interevent times according to theirsizes, implying that only the ordering of interevent timesin T is affected. Each permutation will result in eachrealized event sequence but from exactly the same T .In order to impose power-law distributed burst sizesfor a wide range of ∆ t , we first partition T into severalsubsets, denoted by T l , at different timescales or “levels” l = 0 , , · · · , L : T ≡ { τ i | τ ≤ τ i ≤ ∆ t } ,T l ≡ { τ i | ∆ t l − < τ i ≤ ∆ t l } for l = 1 , · · · , L − , (11) T L ≡ { τ i | τ i > ∆ t L − } , where ∆ t l ≡ τ q s l with constants q , s > l . This partition withsome constants q and s readily determines the numberof bursty trains, denoted by m l , that would be identifiedif ∆ t l were used as the time window, no matter whichpermutation for T is chosen for constructing the eventsequence. It is because each burst size, say b ( l ) , implies b ( l ) − t l and one interevent time larger than ∆ t l . Precisely,we get m l = (cid:12)(cid:12) ∪ Ll (cid:48) = l +1 T l (cid:48) (cid:12)(cid:12) + 1 . (12)Let us now denote the sizes of bursty trains using ∆ t l by B l ≡ { b ( l ) } , with m l = | B l | . Here the burst sizes, b ( l ) s,are to follow a power law with the same exponent β atall levels, for which bursty trains must be hierarchicallyorganized as schematically depicted in Fig. 2. For this,we adopt the power-law burst size distribution only atthe level l = 0: P ∆ t ( b (0) ) = 1 ζ ( β ) b (0) − β for b (0) = 1 , , · · · , (13)where ζ ( · ) denotes the Riemann zeta function. Then m burst sizes are independently drawn from P ∆ t ( b (0) ) inEq. (13) to obtain B . Note that the sum of burst sizesin B must be n + 1. As depicted in Fig. 2, for a given B l ( ( ( (( ( (( ( ((( ( b (0) = b (1) = b (2) = l = FIG. 2. Schematic diagram for the hierarchical organiza-tion of burst trains at different timescales or “levels”, with15 events, denoted by vertical lines, and 14 interevent times.The numbers above the event sequence indicate the levels ofinterevent times, while b ( l ) s are the burst sizes identified using∆ t l . For the definition of ∆ t l , see the main text. with l ≥ B l +1 can be constructed by merging several b ( l ) s to make each b ( l +1) , but under the condition thatthose b ( l +1) s are to be power-law distributed with thesame exponent β in Eq. (13). This condition can be sat-isfied if bigger bursts tend to be merged with bigger ones,and smaller bursts with smaller ones. In other words,bigger (smaller) bursts are followed by bigger (smaller)ones, hence this merging rule can be called the bursty-get-burstier (BGB) method. Precisely, we devise the fol-lowing method: B l is sorted, e.g., in a descending order,then it is sequentially partitioned into m l +1 subsets ofthe (almost) same size. The size of each subset may beeither (cid:98) m l m l +1 (cid:99) or (cid:98) m l m l +1 (cid:99) + 1. The sum of b ( l ) s in eachsubset leads to one b ( l +1) . This procedure is repeateduntil l = L −
1. We numerically confirm that our BGBmethod indeed generates power-law tails in distributionsof b ( l ) s with the same exponent β at all levels, as shownin Fig. 3(a,c,e) for the case with β = 3.We then show how to construct the event sequence x ( t )by permuting interevent times in T , based on informationabout which bursty trains at the level l have been mergedinto which bursty train at the level l + 1 for all possible l s. We begin with the highest level L − L − T L without replacement, denoted by τ ( L ) :( b ( L − , τ ( L )1 , b ( L − , τ ( L )2 , · · · , b ( L − m L − ) . (14)Note that from now on, all the subscripts are dummy in-dexes. Then, each b ( L − is replaced by a sequence madeof b ( L − s, which have been merged together to make the b ( L − , in a random order, alternating with intereventtimes randomly drawn from T L − without replacement.This procedure is repeated for all bursty trains at all lev-els, eventually resulting in a sequence made of exactly n interevent times. From this sequence of interevent times,the event timings are given by t = 0 and t i = (cid:80) ii (cid:48) =1 τ i (cid:48) for i = 1 , · · · , n , then we finally obtain the event sequence x ( t ) by Eq. (8).Once x ( t ) is generated, we first measure distributionsof interevent times and burst sizes for various values of -10 -8 -6 -4 -2 -1 (a) P q ( b ) < b > q b/ q q=1.5361224 β =3 10 -12 -6 -7 P( τ ) -4 -3 -2 -1 -6 -5 -4 -3 -2 -1 (b) α =1.6 A ( t d ) t d β =2.62.72.83.04.0uncor.10 -10 -8 -6 -4 -2 -1 (c) P q ( b ) < b > q b/ q q=1.5361224 β =3 10 -14 -7 -7 P( τ ) -3 -2 -1 -6 -5 -4 -3 -2 -1 (d) α =2 A ( t d ) t d β =2.62.72.83.04.0uncor.10 -10 -8 -6 -4 -2 -1 (e) P q ( b ) < b > q b/ q q=1.5361224 β =3 10 -10 -5 -7 -2 P( τ ) -4 -3 -2 -1 -6 -5 -4 -3 (f) α =2.6 A ( t d ) t d β =2.52.72.93.54.0uncor. FIG. 3. (Color online) Numerical results of the bursty-get-burstier (BGB) model of correlated interevent times for var-ious values of α and β . We show the cases with α = 1 . α = 2 (middle), and 2 . β = 3 for various time windows of∆ t = τ q with q = q s l . (cid:104) b (cid:105) q is the average burst size ofbursty trains identified using ∆ t = τ q . Insets show the cor-responding interevent time distributions. (b), (d), and (f)show the autocorrelation functions for different values of β (colored points), compared to the corresponding uncorrelatedcases (black curves). For each curve, we have generated upto 100 event sequences with n = 5 · , τ = 10 − , q = 1 . s = 2, and L = 4. ∆ t or q = ∆ tτ to confirm our BGB method, and then wecalculate autocorrelation functions to test whether thescaling relation in Eq. (4) holds or not in the presenceof the correlations between interevent times. In Fig. 3,we show the numerical results for various values of α and β , where we have used τ = 10 − , q = 1 . s = 2, and L = 4. We find that P ( τ ) ∼ τ − α and P ∆ t ( b ) ∼ b − β for awide range of ∆ t as expected, see Fig. 5 for more details.We remark that regarding our setting for the power-lawburst size distribution at l = 0 in Eq. (13), one can showthat the value of β is determined for given α , q , and n ,as analyzed in Appendix B. However, by assuming thatit is sufficient to show power-law tails in the burst sizedistributions, we can simulate a wide range of β using γ β α =1.61.92.0 0.2 0.4 0.6 0.8 1 2.5 3 3.5 4(b) γ β FIG. 4. (Color online) The values of γ measured from auto-correlation functions, e.g., in Fig. 3, for various values of α and β , with horizontal dashed lines corresponding to those forthe uncorrelated cases measured in Fig. 1(b). our BGB method.Then, for each α , autocorrelation functions A ( t d ) forvarious values of β are compared to that for the uncorre-lated case, e.g., in Fig. 3(b,d,f). The estimated values of γ for various values of α and β are presented in Fig. 4.When α ≤
2, it is numerically found that the autocorre-lation functions for β (cid:46) α and γ in Eq. (4). Precisely, the smaller β leads to thelarger γ , implying that the stronger correlations betweeninterevent times may induce the faster decaying of auto-correlation. The deviation observed for β (cid:46) b diverges for β < α >
2, the estimated γ deviates significantly from that for the uncorrelated casefor the almost entire range of β , although γ seems toapproach the uncorrelated case as β increases. Interest-ingly, the estimated values of γ show an increasing andthen decreasing behavior as β increases. The reason forsuch different behaviors of γ for α ≤ α > γ as afunction of α in the uncorrelated cases. In order to un-derstand this difference, more rigorous studies are neededin a future. III. CONCLUSION
The effects of correlations between interevent times,often called correlated bursts, on the temporal correla-tions characterized by power-law decaying autocorrela-tion functions, are far from being fully understood. Inorder to study these effects systematically, we have de-vised the bursty-get-burstier (BGB) model, where power-law distributions of interevent times and burst sizes areinputs rather than outputs of the model. With ourmodel, one can tune the degree of correlations betweeninterevent times, while keeping the same interevent timestatistics. Then the established scaling relation betweenpower-law exponent α for interevent time distributionsand decaying exponent γ for autocorrelation functionscan be numerically tested, especially, whether the scaling relation can be violated due to the correlations betweeninterevent times.As a baseline, we numerically study the case of un-correlated interevent times. We find that the depen-dence of γ on α is comparable to the theoretical expec-tation in Eq. (4) for the range of α far from 2. It isbecause for α ≈
2, the logarithmic corrections becomeeffective. Next, using our BGB model, we generate theevent sequences showing power-law distributions of bothinterevent times and burst sizes for a wide range of thetime window. By measuring the autocorrelation func-tions of those event sequences and then by estimatingthe decaying exponents of them, we find that the corre-lations between interevent times can violate the scalingrelation established for the uncorrelated case, but in dif-ferent ways depending on the range of α .Our BGB model for correlated bursts turns out to beuseful for imposing power-law distributed burst sizes atdifferent timescales simultaneously, which is however nottrivial to implement. It is because it requires somewhatdeliberate hierarchical organization of burst sizes at dif-ferent timescales, namely the bursty-get-burstier mech-anism. There can be other ways of implementing suchhierarchical organization. We can get some hints fromthis hierarchical organization for the origin of correlatedbursts. For example, if human communication patternscan be described in terms of correlated bursts, humanindividuals might have organized their communicationactivities in a hierarchical way either consciously or un-consciously: Bigger bursts tend to be followed by biggerones, while smaller bursts follow smaller ones, at all rel-evant timescales of human dynamics. ACKNOWLEDGMENTS
The author thanks Mikko Kivel¨a, J´anos Kert´esz, andKimmo Kaski for fruitful discussions, and he acknowl-edges financial support by Basic Science Research Pro-gram through the National Research Foundation of Ko-rea (NRF) grant funded by the Ministry of Education(2015R1D1A1A01058958).
Appendix A: Asymptotic result of the convolution ofthe interevent time distribution
In Eq. (10), P (cid:63)k ( t d ) denotes the k th order convolutionof the interevent time distribution P ( τ ), i.e., P (cid:63) ( t d ) = P ( τ = t d ) and P (cid:63)k ( t d ) = (cid:82) t d P ( τ ) P (cid:63)k − ( t d − τ ) dτ for k >
1. In other words, one can write as follows [32]: P (cid:63)k ( t d ) = k (cid:89) i =1 (cid:90) ∞ dτ i P ( τ i ) · δ (cid:32) t d − k (cid:88) i =1 τ i (cid:33) . (A1) -12 -10 -8 -6 -4 -2 -1 (a) α =1.6 P q ( b ) < b > q b/ q q=1.5361224 β =2.4 10 -12 -10 -8 -6 -4 -2 -1 (b) P q ( b ) < b > q b/ q q=1.5361224 β =2.6 10 -10 -8 -6 -4 -2 -1 (c) P q ( b ) < b > q b/ q q=1.5361224 β =3 10 -10 -8 -6 -4 -2 -1 (d) P q ( b ) < b > q b/ q q=1.5361224 β =410 -12 -10 -8 -6 -4 -2 -1 (e) α =2 P q ( b ) < b > q b/ q q=1.5361224 β =2.4 10 -12 -10 -8 -6 -4 -2 -1 (f) P q ( b ) < b > q b/ q q=1.5361224 β =2.6 10 -10 -8 -6 -4 -2 -1 (g) P q ( b ) < b > q b/ q q=1.5361224 β =3 10 -10 -8 -6 -4 -2 -1 (h) P q ( b ) < b > q b/ q q=1.5361224 β =410 -12 -10 -8 -6 -4 -2 -1 (i) α =2.6 P q ( b ) < b > q b/ q q=1.5361224 β =2.4 10 -10 -8 -6 -4 -2 -1 (j) P q ( b ) < b > q b/ q q=1.5361224 β =2.6 10 -10 -8 -6 -4 -2 -1 (k) P q ( b ) < b > q b/ q q=1.5361224 β =3 10 -10 -8 -6 -4 -2 -1 (l) P q ( b ) < b > q b/ q q=1.5361224 β =4 FIG. 5. (Color online) Burst size distributions of the bursty-get-burstier model of correlated interevent times for α = 1 . . β = 2 .
4, 2 .
6, 3, and 4 (from left to right), and for various values of ∆ t = τ q with q = q s l . (cid:104) b (cid:105) q denotes the average burst size of bursty trains identified using ∆ t = τ q . For each case, we have generated 20 eventsequences with n = 5 · , τ = 10 − , q = 1 . s = 2, and L = 4. Here we analyze P (cid:63)k ( t d ) for t d (cid:29) τ . By taking a Laplacetransform of P (cid:63)k ( t d ) with respect to t d , one gets (cid:103) P (cid:63)k ( s ) = ˜ P ( s ) k , (A2)where ˜ P ( s ) denotes the Laplace transform of P ( τ ) inEq. (9) and is given by˜ P ( s ) = ( α − τ s ) α − Γ(1 − α, τ s ) . (A3)The incomplete Gamma function is expanded in theasymptotic limit of s → P ( s ) ≈ − α )( α − τ s ) α − − α − α − τ s + · · · . (A4)If 1 < α <
2, since the term of the order of s α − domi-nates that of s , we can ignore the higher order terms toget ˜ P ( s ) k ≈ k Γ(1 − α )( α − τ s ) α − (A5)= 1 + Γ(1 − α )( α − τ k s ) α − , (A6)where we have defined τ k by τ α − k ≡ kτ α − . That is, the k th convolution can interpreted as the replacement of τ by τ k . We finally get for t d (cid:29) τ P (cid:63)k ( t d ) ≈ ( α − τ α − k t − αd (A7)= k ( α − τ α − t − αd , (A8)which turns out to increase as α varies from 1. If α > s in Eq. (A4) becomesdominant, we can similarly obtain˜ P ( s ) k ≈ − k α − α − τ s, (A9) P (cid:63)k ( t d ) ≈ k α − ( α − τ α − t − αd . (A10)This P (cid:63)k ( t d ) must be decreasing according to α as long as kτ t d <
1. In sum, the convolution of the interevent timedistribution is increasing for small α and decreasing forlarge α , implying the non-monotonic behavior of P (cid:63)k ( t d )according to α . Appendix B: Exact scaling relation between α and β Here we show that in principle, α and β cannot beindependent of each other. As mentioned in the maintext, the sum of burst sizes in B must be equal to thetotal number of events, i.e., n + 1, implying that m (cid:104) b (0) (cid:105) = n + 1 , (B1)where (cid:104) b (0) (cid:105) denotes the average burst size. Note thatthis equation holds for the arbitrary choice of P ( τ ) and P ∆ t ( b (0) ). In our setting with Eqs. (9) and (13), since m = n (cid:90) ∞ ∆ t P ( τ ) dτ = nq − α , (B2) (cid:104) b (0) (cid:105) = ∞ (cid:88) b (0) =1 b (0) P ∆ t ( b (0) ) = ζ ( β − ζ ( β ) , (B3) we obtain the relation between α and β from Eq. (B1) asfollows: ζ ( β − ζ ( β ) = n + 1 n q α − . (B4)This result can be seen as another scaling relation be-tween α and β for given q and n . For example, when q = 1 . n (cid:29)
1, one gets β ≈ .
81 for α = 2, and β ≈ .
51 for α = 2 .
6, respectively. We remark that thisrelation is based on the assumption that the burst sizedistribution follows a clear power law for the entire rangeof b (0) . On the other hand, in practice, we can simulatea much wider range of β by considering the distributionsshowing power laws only in their tails. [1] M. S. Wheatland, P. A. Sturrock, and J. M. McTiernan,The Astrophysical Journal , 448 (1998).[2] A. Corral, Physical Review Letters , 108501 (2004).[3] L. de Arcangelis, C. Godano, E. Lippiello, andM. Nicodemi, Physical Review Letters , 051102 (2006).[4] T. Kemuriyama, H. Ohta, Y. Sato, S. Maruyama,M. Tandai-Hiruma, K. Kato, and Y. Nishida, BioSys-tems , 144 (2010).[5] A.-L. Barab´asi, Nature , 207 (2005).[6] T. Nakamura, K. Kiyono, K. Yoshiuchi, R. Nakahara,Z. R. Struzik, and Y. Yamamoto, Physical Review Let-ters , 138103 (2007).[7] P. Bak, C. Tang, and K. Wiesenfeld, Physical ReviewLetters , 381 (1987).[8] M. B. Weissman, Reviews of Modern Physics , 537(1988).[9] L. Ward and P. Greenwood, Scholarpedia , 1537 (2007).[10] M. Karsai, K. Kaski, A.-L. Barab´asi, and J. Kert´esz,Scientific Reports , 397 (2012).[11] P. Panzarasa and M. Bonaventura, Physical Review E , 062821 (2015).[12] A. Vazquez, B. R´acz, A. Luk´acs, and A. L. Barab´asi,Physical Review Letters , 158702 (2007).[13] M. Karsai, M. Kivel¨a, R. K. Pan, K. Kaski, J. Kert´esz,A.-L. Barab´asi, and J. Saram¨aki, Physical Review E ,025102 (2011).[14] G. Miritello, E. Moro, and R. Lara, Physical Review E , 045102 (2011).[15] L. E. C. Rocha, F. Liljeros, and P. Holme, PLoS Com-putational Biology , e1001109 (2011).[16] H.-H. Jo, J. I. Perotti, K. Kaski, and J. Kert´esz, PhysicalReview X , 011041 (2014).[17] J.-C. Delvenne, R. Lambiotte, and L. E. C. Rocha, Na-ture Communications , 7366 (2015).[18] S. B. Lowen and M. C. Teich, Physical Review E , 992(1993). [19] S. Vajna, B. T´oth, and J. Kert´esz, New Journal ofPhysics , 103023 (2013).[20] S. Abe and N. Suzuki, Physica A: Statistical Mechanicsand its Applications , 1917 (2009).[21] M. Karsai, K. Kaski, and J. Kert´esz, PLoS ONE ,e40612 (2012).[22] W. Wang, N. Yuan, L. Pan, P. Jiao, W. Dai, G. Xue,and D. Liu, Physica A: Statistical Mechanics and its Ap-plications , 846 (2015).[23] Z.-Q. Jiang, W.-J. Xie, M.-X. Li, W.-X. Zhou, andD. Sornette, Journal of Statistical Mechanics: Theoryand Experiment , 073210 (2016).[24] H.-H. Jo, J. I. Perotti, K. Kaski, and J. Kert´esz, PhysicalReview E , 022814 (2015).[25] P. Allegrini, D. Menicucci, R. Bedini, L. Fronzoni,A. Gemignani, P. Grigolini, B. J. West, and P. Para-disi, Physical Review E , 061914 (2009).[26] J. W. Kantelhardt, E. Koscielny-Bunde, H. H. A. Rego,S. Havlin, and A. Bunde, Physica A: Statistical Mechan-ics and its Applications , 441 (2001).[27] D. Rybski, S. V. Buldyrev, S. Havlin, F. Liljeros, andH. A. Makse, Proceedings of the National Academy ofSciences , 12640 (2009).[28] D. Rybski, S. V. Buldyrev, S. Havlin, F. Liljeros, andH. A. Makse, Scientific Reports , 560 (2012).[29] J. I. Perotti, H.-H. Jo, P. Holme, and J. Saram¨aki, “Tem-poral network sparsity and the slowing down of spread-ing,” (2014), arXiv:1411.5553.[30] E.-K. Kim and H.-H. Jo, Physical Review E , 032311(2016).[31] In order to avoid too large interevent times, especiallyfor small α , we introduce the upper bound of intereventtimes for numerical simulations. This upper bound is setto be 1, hence being believed to have a negligible effecton the results.[32] H.-H. Jo, R. K. Pan, J. I. Perotti, and K. Kaski, PhysicalReview E87