The spreading of computer viruses on time-varying networks
TThe spreading of computer viruses on time-varying networks
Terry Brett, George Loukas, Yamir Moreno,
2, 3 and Nicola Perra
1, 3, ∗ University of Greenwich, Old Royal Naval College, London, UK Institute for Biocomputation and Physics of Complex Systems (BIFI), University of Zaragoza, Zaragoza, Spain ISI Foundation, Turin, Italy (Dated: January 23, 2019)Social networks are the prime channel for the spreading of computer viruses. Yet the study of their propa-gation neglects the temporal nature of social interactions and the heterogeneity of users’ susceptibility. Here,we introduce a theoretical framework that captures both properties. We study two realistic types of virusespropagating on temporal networks featuring Q categories of susceptibility and derive analytically the invasionthreshold. We found that the temporal coupling of categories might increase the fragility of the system to cyberthreats. Our results show that networks’ dynamics and their interplay with users features are crucial for thespreading of computer viruses.
Alongside clear societal and economic benefits, moderntechnology exposes us to serious challenges. In particular, thespreading of malicious content online, often based on inge-nious deception strategies, is one of the most pressing becauseit poses serious threats to our privacy, finances, and safety [1].Victims of a typical social engineering attack [2] may receivea message containing a malicious link or file, appearing tooriginate from a friend or other trusted entity. If opened, itmay compromise the computer, access personal information,and spread the virus further unbeknownst to the victim. Re-cent research has shown how the susceptibility of individualsto such attacks is not homogenous and depends on several fea-tures such as age, prior training, computer proficiency, famil-iarity with social network platforms, among others [3–5]. Fur-thermore, the properties of real networks are known to facili-tate the propagation of such processes [6–15]. In particular,the heterogeneity in contact patterns makes socio-technicalsystems quite fragile to biological and digital threats.The study of these phenomena has largely neglected thecomplex temporal nature of online contact patterns in favorof static and time-aggregated approaches [16, 17]. These ap-proximations might be fitting. Indeed, in the past, computerviruses would spread mainly via email networks, targeting theaddress books of victims, which contain contacts lists [18].However, not many people create such lists any more and ac-cess to them is restricted [7]. In the context of social or bio-logical contagions, neglecting the temporal nature of the net-works where the processes unfold has been shown to inducemisrepresentations of their spreading potential. In fact, theorder and concurrency of connections is key [19–43]. To thebest of our knowledge, beside some early work on the spread-ing of viruses via Bluetooth among mobile phones [44], thestudy of the propagation of cyber threats considering the tem-poral nature of social interactions is still missing. Further-more, with few exceptions [45], the literature devoted to thestudy of computer viruses unfolding on networks typically ne-glects that the susceptibility of online users is not homoge-nous. Conversely, the literature that studies the susceptibilityof users to cyber threats traditionally focuses on single usersneglecting their connections.To tackle these limitations, here we introduce a theoretical framework to study the spreading of computer viruses, basedon social engineering deception strategies, on time-varyingnetworks. We model users’ interactions using a time-varyingnetwork model and consider two types of viruses. The firstmimics threats that can propagate only via connections acti-vated at each time step. The second, on the contrary, considersviruses able to access also information about past connections.We investigate the impact of different classes of susceptibilityconsidering that they might also influence the link formationprocess. In all cases, we analytically derive the conditions reg-ulating the spreading of the virus. Interestingly, these are de-fined by the interplay between the features of the cyber threats,the categories of susceptibility and their time-varying connec-tivity. Furthermore, in some scenarios, the coupling betweencategories creates a complex phenomenology that favors thespreading of the virus. These results have the potential to ini-tiate future efforts aimed at describing more realistically thespreading of computer viruses on online social networks.We consider a population of N online users which ex-change messages in a time-varying network. Nodes are as-signed to one of Q categories describing their susceptibility tocyber threats measured in terms of their gullibility and timeneeded to recover from successful attacks. Since suscepti-bility is linked to demographic features, we consider that themembership to a category might influence the link creationprocess. In fact, homophily is a strong social mechanismknown to affect the structure and organization of ties [46].We model the contact patterns between users with a general-ization of the activity-driven framework [21, 47–49]. Here,nodes feature an activity a describing their propensity to initi-ate communications. Activities are extracted from a distribu-tion F ( a ) which, as observations in real systems have shown,is typically heterogenous [21, 22, 48, 50]. We select power-law distributions F ( a ) ∼ a − α with a ∈ [ (cid:15), to avoid diver-gences. At each time step nodes are active with probability a ∆ t . Active nodes select m others and create directed (out-going) links which mimic messages.In the simplest version of activity-driven networks the se-lection is random and memoryless [21]. Here, we proposea variation: with probability p each target is selected, at ran-dom, among the group of nodes in the same category, and with a r X i v : . [ phy s i c s . s o c - ph ] J a n probability − p among the nodes in any other category. Inother words, p tunes the homophily level in the network withrespect to susceptibility to cyber threats. At time t + ∆ t alledges are deleted and the process starts from the beginning.Unless specified otherwise, links have a duration ∆ t . With-out loss of generality we set ∆ t = 1 . The model is clearlya simplification of real interactions. However, it offers sim-ple, yet non trivial, settings to study the effects of temporalconnectivity patterns on contagion processes unfolding at acomparable time-scale with respect to the evolution of con-nections [20, 21, 47, 51].We describe the propagation of a computer virus adoptingthe prototypical SIS model [13, 52]. At each time step t thevirus, unbeknownst to the victims, sends a message, with ma-licious content, to all the nodes genuinely contacted at t (virustype 1) or within t − τ time-step (virus type 2). The focus isnot defining the optimal set of nodes to maximize/minimizethe damage. Thus, we select randomly a small percentage( . ) of nodes as initial seeds. In these settings, suscepti-ble nodes of class x ∈ [1 , . . . , Q ] , that receive a maliciousmessage, become infectious with probability λ x which de-fines their gullibility. They recover and become susceptibleagain with rate µ x . Assuming that nodes with the same valueof activity and in the same category are statistically equiva-lent, we group nodes according to the two features. At eachtime step, we call S xa and I xa the number of nodes suscepti-ble and infected in activity class a and category x . Clearly (cid:82) daS xa = S x , (cid:82) daI xa = I x , (cid:80) x S x = S , and (cid:80) x I x = I .Furthermore, N xa describes the number of nodes of activity a in category x , thus (cid:82) daN xa = N x and (cid:80) x N x = N . Inthese settings, we can represent the variation of the number ofinfected nodes of activity a in category x as: d t I xa = − µI xa + λ x mS xa × p (cid:90) da (cid:48) a (cid:48) I xa (cid:48) N x + (1 − p ) (cid:88) y (cid:54) = x (cid:90) da (cid:48) a (cid:48) I ya (cid:48) N − N y . (1)The first term on the right hand side accounts for the recov-ery process. The second and third terms capture susceptiblenodes that receive messages from active and infected verticesin the same (second) or different (third) category, and get in-fected as a result. With respect to the typical biological con-tagion process, here transmission is asymmetric. Only nodesreceiving a message from an infected person might be exposedto the virus. Thus, not only the order of connections, butalso their direction is a crucial ingredient for the spreading.Since the links are created randomly, each node is selectedwith a probability pm/N x by nodes in the same category or (1 − p ) m/ ( N − N y ) by nodes in other categories. The totalnumber of nodes is constant thus S xa = N xa − I xa and at theearly stages of the spreading we can assume that the numberof infected nodes is very small: S xa ∼ N xa . By integrating across all activities Eq. 3 we get: d t I x = − µ x I x + λ x m pθ x + (1 − p ) N x (cid:88) y (cid:54) = x θ y / ( N − N y ) , where we define θ x = (cid:82) daaI xa . By multiplying both sides ofEq. 3 for a and integrating across all the activities we obtain d t θ x = − µ x θ x + mλ x (cid:104) a (cid:105) x pθ x + (1 − p ) N x (cid:88) y (cid:54) = x θ y / ( N − N y ) . The virus is able to spread, if and only if the largest eigenvalueof the Jacobian matrix of the system of differential equationsin I x and θ x is larger than zero [21]. As shown in details inthe Supplementary Material (SM) this implies: R = p (cid:80) x β x + Ξ (cid:80) x µ x > , (2)where R is the basic reproductive number defined as the av-erage number of infected nodes generated, in a fully suscepti-ble population, by an infected individual [52], β x = mλ x (cid:104) a (cid:105) x and Ξ is a function of the interplay between the average acti-vation, infection and recovery rate of each category as well asof the mixing between categories.To understand the dynamics, let us consider a particularcase in which the system is characterized by only two cat-egories. Furthermore, let us consider, as first scenario, thatall nodes have the same recovery rate. In these settings wehave Ξ = p ( β + β ) + 4 β β (1 − p ) . The conditionfor the spreading, even with only two classes, is a non linearfunction of the average activity of each category, the infec-tion probabilities per contact and the homophily. In the limit p = 0 , nodes in a category connects only with vertices inthe other and the expression reduces to R = √ β β µ . In thelimit p = 1 instead, interactions are only between nodes inthe same category. The system is effectively split in two dis-connected networks and there are two independent conditions R x = β x /µ . For a general p , these two values confine R : min x R x ≤ R ( p ) ≤ max x R x . In fact, any value of p < will reduce the spreading power of the category characterizedwith the largest R x as some connections will be establishedwith nodes where the virus finds it harder to spread (see SMfor the proof).In Fig. 1-A-C, we compare analytical predictions with nu-merical simulations. We set λ = 0 . and use Eq. 2 to esti-mate the critical value of λ for which R ≡ . On the y -axiswe plot the lifetime of the process defined as the time that thevirus needs either to die out or to reach a fraction Y of the pop-ulation [53]. The lifetime acts as the susceptibility of a secondorder phase transition and allows a precise numerical estima-tion of the threshold of SIS processes [53]. In panels A-B weconsider a scenario in which nodes are assigned randomly toone of the two categories. Thus the average activity in the two A) B) C)D) E) F) λ FIG. 1. Lifetime of the SIS process (A-C) and contour plot of R (D-F). In A-B-D-E nodes are randomly assigned to two categories,in C-F instead in decreasing order of activity. We set p = 0 . (A-D), p = 0 . (B-C-E-F). In A-C we fix N = 2 × , m = 4 , α = − . , µ = µ = 10 − , λ = 0 . , Y = 0 . , and . of randominitial seeds. We plot the median and confidence intervals in simulations per point. The solid lines come from Eq. 2, and thedashed lines are the analytical threshold in case of a single category.In the contour plot we set µ = µ = 10 − . is the same and set p = 0 . and p = 0 . respectively. Theanalytical value of the threshold (vertical solid line) perfectlymatches the numerical estimation. For p = 0 . the thresholdis smaller than for p = 0 . and closer to the threshold of asystem with a single category (dashed lines). For smaller val-ues of homophily, instead, the critical conditions are driven bythe interplay between the activation rates and gullibility of thetwo categories. Panels D-E show the analytical value of R as a function of λ and λ for the two values of p . The greyregions are sub-critical, i.e., the virus is not able to spread.Since the average activity in the two categories is the same,the two plots are symmetric. Interestingly, the active region(where the virus is able to spread) is larger for large values of p . This is due to the fact that in these settings the virus willspread if above the threshold in at least one category inde-pendently of the other. In the opposite limit, on the contrary,the two categories get intertwined and a small value of the in-fection probability in one category should be associated to aprogressively large value in the other.In panels C-F we consider that the first category containsa fraction g of nodes selected in decreasing order of activ-ity. Thus, this category contains the gN most active nodes,while the other the (1 − g ) N least active (see SM). To com-pare with panel B, we set g = 0 . and p = 0 . . First, the an-alytical threshold nicely matches the numerical simulations.Second, although the other parameters are the same used inpanel B, the critical value of the gullibility of the first class issmaller. Thus, correlations between activity and gullibility fa-cilitate the spreading. This is confirmed in panel F where theactive phase space features a region in which the spreading iscompletely dominated by the category of most active nodes.Overall, all the plots show the importance of distinguishingnodes according to their gullibility. Indeed, neglecting thepresence of different classes of users might induce a strongmisrepresentation of the virus propagation (dashed lines).Let us next consider a second scenario where categories dif-ferentiate also for the time needed to recover from a successful attack. For two categories, we can write Ξ = ( µ − µ ) + p ( β + β ) + 2 p ( µ − µ )( β − β ) + 4 β β (1 − p ) .Interestingly, we have the same terms that appeared in thefirst scenario, plus two that feature the difference betweenthe recovery rates and β s of the two categories. Thus R isa function of the interplay between the activities, gullibili-ties and recovery rates. In the limit p = 0 , each categoryonly connects with nodes in the other, the two groups are cou-pled and the threshold reads R = √ ( µ − µ ) +4 β β µ + µ . In thelimit p = 1 instead, the two categories are completely de-coupled and the threshold becomes, as before, R = β x /µ x .As shown in Fig. 2-D-H, for a general value of p the repro-ductive number is not bounded, as before, by the values of R x computed in the two classes separately (see SM). In Fig. 2-D, we assign nodes randomly to each category, fix β x and µ x and compute R as a function of p . In the shaded area min x R x ≤ R ( p ) ≤ max x R x . Interestingly, after a p ∗ (ver-tical dashed line), which as shown in the SM can be computedanalytically, we enter in a regime where R ( p ) > max x R x .Thus, only specific values of the coupling between categoriesmight induce the virus to spread faster in the combined sys-tem than in each single category in isolation. However, thisnon linear effect is found only in a small fraction of the phasespace see Fig. 2-H. The necessary, but not sufficient condi-tion, is that two categories differentiate both for gullibility andrecovery rates in such a way that one is more gullible and re-covers faster than the other. In this regime, the right mixingbetween the two might create a feedback loop that makes thesystem more fragile.Fig. 2-A-C shows a good match between the analytical(solid vertical lines) and numerical thresholds in case of nodesare assigned at random (A-B) or in decreasing order of activ-ity (C) to the two categories. We fix two different recoveryrates, λ , and use λ as order parameter. Panels A-B dif-fer in the value of the homophily p . We set p = 0 . in A,while p = 0 . in B-C. The presence of a category of nodescharacterized by a smaller value of recovery rate pushes thethreshold to smaller values with respect to the first scenarios(Fig. 1). As before, the value of the threshold estimated con-sidering only a single category, characterized by the averagerecovery rate of the two, (dashed lines) leads to a misrepre-sentation of the spreading power of the virus, especially forsmaller values of homophily (see panel B).The effect of p on the critical value of λ is similar to thefirst scenario. In fact, even when categories differentiate bythe recovery rates, high values of homophily push the criticalpoint to smaller values. However, here the difference betweenthe two is less significant than in Fig. 1. In Fig. 2-E-F, weshow the analytical value of R as function of µ and µ . In-terestingly, the sub-critical region, for p = 0 . , is smaller thanfor p = 0 . . This is in contrast to what was observed in thecorresponding plots for the first scenario and highlights onceagain the complex phenomenology introduced by the inter-play of different recovery rates. In Fig. 2-C-G we investigatethe effect of correlations. In case that the most active nodes A) B) C) D)E) F) G) H) μ FIG. 2. Lifetime of the process (A-C), R ( µ , µ ) (E-G), R ( p ) (D), and p ∗ ( µ , µ ) (H). In A-B-E-F nodes are randomly assignedto two categories, in C-G instead in decreasing order of activity. Weset p = 0 . (A-E), p = 0 . (B-C-F-G). In panels A-C we set N =2 × , m = 4 , α = − . , µ = 10 − , µ = 5 × − , λ = 0 . , Y = 0 . , and . randomly selected seeds. We plot the medianand confidence intervals in simulations per point. The solidlines come from Eq. 2. The dashed lines are the analytical thresholdin case of a single category of recovery rate characterized by theaverage value of the recovery rates. In the contour plot we set λ =0 . and λ = 0 . . In D the shaded area describe the region where min x β x /µ x ≤ R ≤ max x β x /µ x . The dashed vertical line thethe analytical value of p above which R > max x β x /µ x . In H weplot p ∗ as function of µ and µ . In D-H we set λ = 0 . and usethe same parameters of the other plots. are able to recover quickly from the attack, the virus is ableto spread only if the gullibility of such users is higher than inthe corresponding case without correlations (panel B). This isconfirmed in panel G, where we see that correlations betweenrecovery rates significantly change the active region.Finally, we turn our attention to a second type of virusable to access also past contacts of infected users within atime window τ . As before, the virus propagates via activeinfected nodes, but at each time t active users might infecttheir contacts in a time-window ( t − τ, t ] . Within a mean-fieldapproximation, we can adopt the same equations describedabove and change the probability that a node in each activ-ity class receives a message by active and infected nodes. Inthis case, the out-degree of each active node is not m , but afunction of τ : k out ( a ) = m (cid:2) a + ( τ − a (cid:3) (see SM). Tograsp the derivation, consider the simplest scenario in which τ = 2 . In this case, active nodes might have either m or m contacts in two time steps. The first class describesnodes that are active at time t but were not active at time t − ; whereas the second, nodes that were active in bothtime steps. Thus the out-degree of these nodes, on average, is k out ( a ) = ma (1 − a ) + 2 ma . As shown in the SM, the con-dition for the spreading has the same structure of Eq. 2 where,however, the value of β s are changed with the following trans-formation m → m (cid:2) (cid:104) a (cid:105) + ( τ − (cid:104) a (cid:105) (cid:3) . Thus, the larger thevisibility of past connections, from the virus point of view, thelarger R . Intuitively this is due to the fact that the virus, forlarge values of τ , is able to access more contacts, which resultsin a larger spreading potential. This observation nicely showshow neglecting the temporal nature of connectivity patterns infavor of static (or time integrated) approximations might leadto a poor description of the propagation of viruses that do not λ C λ C λ C A) B) C)
FIG. 3. Lifetime of the SIS process for τ = 2 , , (A,B,C) for twocategories to which nodes are assigned randomly. Simulations aredone setting N = 2 × , m = 4 , α = − . , Y = 0 . , µ = 10 − , λ = 0 . , p = 0 . , and . random initial seeds. We plot themedian and confidence intervals in simulations per point. have access to contacts lists or past connections. In Fig. 3 weshow the comparison between analytical (solid lines) and nu-merical values of the threshold for different values of τ . Toisolate the effect of τ we considered two categories, a singlerecovery rate, and set p = 0 . . The analytical value is a goodapproximation only for small values of τ . The mean-field ap-proximation becomes less accurate as more connections frompast time-steps are kept in memory. Thus, the analytical esti-mation provides only a lower bound, which together with thesolution for τ = 1 (dashed lines) − that constitutes an upperbound − , marks the region where spreading is possible (red re-gions). In other words, for a general value of τ , the thresholdwill be lower than the analytical value computed for τ = 1 ,and larger than the corresponding value computed at τ .Overall our results highlight how the spreading of computerviruses based on social engineering is critically affected bythe temporal nature of our interactions and different suscep-tibilities to cyber threats. Our findings show that networks’dynamics and their interplay with the characteristics of usershave to be considered in order to avoid misrepresentation ofthe spreading power of computer viruses in social networks.We have also quantified the extent to which the previous mis-match is important for three plausible scenarios. We, how-ever, note that we have studied a simple network model thatneglects a range of properties of real social networks such asthe presence of weak and strong ties, high order correlations,and community structures. The study of the impact of thesefeatures on the unfolding of computer viruses calls for addi-tional research.This material is based upon work supported by, or in partby, the U. S. Army Research Laboratory and the U. S. ArmyResearch Office under contract/grant number W911NF-18-1-0376. Y. M. acknowledges support from the Governmentof Arag´on, Spain through grant E36-17R (FENOL) and byMINECO and FEDER funds (grant FIS2017-87519-P). Theauthors thanks Andrea Baronchelli and Michele Starnini foruseful discussions. ∗ [email protected][1] I. Kayes and A. Iamnitchi, Online Social Networks and Media , 1 (2017). [2] R. Heartfield and G. Loukas, ACM Computing Surveys(CSUR) , 37 (2016).[3] R. Heartfield and G. Loukas, Computers & Security , 101(2018).[4] R. Heartfield, G. Loukas, and D. Gan, in IEEE 15th Interna-tional Conference on Software Engineering Research, Manage-ment and Applications (SERA) (IEEE, 2017) pp. 371–378.[5] R. Heartfield, G. Loukas, and D. Gan, IEEE Access , 6910(2016).[6] A. L. Lloyd and R. M. May, Science , 1316 (2001).[7] J. Balthrop, S. Forrest, M. E. Newman, and M. M. Williamson,Science , 527 (2004).[8] R. Pastor-Satorras and A. Vespignani, Phys. Rev, Lett. , 3200(2001).[9] Y. Moreno and A. V´azquez, Eur. Phys. J. , 265 (2003).[10] M. E. J. Newman, Phys. Rev. E , 016128 (2002).[11] M. Newman, Networks. An Introduction (Oxford UnivesityPress, 2010).[12] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, andA. Vespignani, Reviews of Modern Physics , 925 (2015).[13] A. Barrat, M. Barth´elemy, and A. Vespignani, DynamicalProcesses on Complex Networks (Cambridge Univesity Press,2008).[14] L.-X. Yang, X. Yang, J. Liu, Q. Zhu, and C. Gan, AppliedMathematics and Computation , 8705 (2013).[15] L.-X. Yang and X. Yang, Physica A: Statistical Mechanics andIts Applications , 173 (2014).[16] P. Holme, The European Physical Journal B , 1 (2015).[17] P. Holme and J. Saram¨aki, Physics Reports , 97 (2012).[18] M. E. Newman, S. Forrest, and J. Balthrop, Physical Review E , 035101 (2002).[19] A. Barrat and C. Cattuto, in Social Phenomena (Springer Inter-national Publishing, 2015) pp. 37–57.[20] N. Perra, A. Baronchelli, D. Mocanu, B. Gonc¸alves, R. Pastor-Satorras, and A. Vespignani, Physical Review Letter ,238701 (2012).[21] N. Perra, B. Gonc¸alves, R. Pastor-Satorras, and A. Vespignani,Scientific Reports , 469 (2012).[22] B. Ribeiro, N. Perra, and A. Baronchelli, Scientific Reports ,3006 (2013).[23] S. Liu, N. Perra, M. Karsai, and A. Vespignani, Physical Re-view Letters , 118702 (2014).[24] S.-Y. Liu, A. Baronchelli, and N. Perra, Physical Review E ,032805 (2013).[25] G. Ren and X. Wang, Chaos: An Interdisciplinary Journal ofNonlinear Science , 023116 (2014).[26] M. Starnini, A. Machens, C. Cattuto, A. Barrat, and R. Pastor-Satorras, Journal of Theoretical Biology , 89 (2013).[27] M. Starnini, A. Baronchelli, A. Barrat, and R. Pastor-Satorras,Physical Review E , 056115 (2012).[28] E. Valdano, L. Ferreri, C. Poletto, and V. Colizza, PhysicalReview X , 021005 (2015).[29] I. Scholtes, N. Wider, R. Pfitzner, A. Garas, C. Tessone, andF. Schweitzer, Nature Communications , 5024 (2014).[30] M. J. Williams and M. Musolesi, Royal Society Open Science , 160196 (2016).[31] L. E. Rocha and N. Masuda, New Journal of Physics , 063023(2014).[32] T. Takaguchi, N. Sato, K. Yano, and N. Masuda, New Journalof Physics , 093003 (2012).[33] L. E. Rocha and V. D. Blondel, PLoS computational biology ,e1002974 (2013).[34] G. Ghoshal and P. Holme, Physica A: Statistical Mechanics andits Applications , 603 (2006). [35] K. Sun, A. Baronchelli, and N. Perra, The European PhysicalJournal B , 1 (2015).[36] D. Mistry, Q. Zhang, N. Perra, and A. Baronchelli, PhysicalReview E , 042805 (2015).[37] R. Pfitzner, I. Scholtes, A. Garas, C. Tessone, andF. Schweitzer, Physical Review Letter , 19 (2013).[38] T. Takaguchi, N. Sato, K. Yano, and N. Masuda, New Journalof Physics , 093003 (2012).[39] T. Takaguchi, N. Masuda, and P. Holme, PloS one , e68629(2013).[40] P. Holme and F. Liljeros, Scientific Reports , 4999 (2014).[41] P. Holme and N. Masuda, PloS one , e0120567 (2015).[42] Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Man-fredi, M. Perc, N. Perra, M. Salath´e, and D. Zhao, PhysicsReports , 1 (2016).[43] B. Gonc¸alves and N. Perra, Social phenomena: From data anal-ysis to models (Springer, 2015).[44] P. Wang, M. C. Gonz´alez, C. A. Hidalgo, and A.-L. Barab´asi,Science , 1071 (2009).[45] S. Peng, G. Wang, Y. Zhou, C. Wan, C. Wang, and S. Yu, IEEETransactions on Dependable and Secure Computing (2017).[46] M. McPherson, L. Smith-Lovin, and J. M. Cook, Annual re-view of sociology , 415 (2001).[47] M. Karsai, N. Perra, and A. Vespignani, Scientific Reports ,4001 (2014).[48] E. Ubaldi, N. Perra, M. Karsai, A. Vezzani, R. Burioni, andA. Vespignani, Scientific Reports , 35724 (2016).[49] M. Tizzani, S. Lenti, E. Ubaldi, A. Vezzani, C. Castellano, andR. Burioni, Physical Review E , 062315 (2018).[50] M. Tomasello, N. Perra, C. Tessone, M. Karsai, andF. Schweitzer, Scientific Reports , 5679 (2014).[51] S. Liu, N. Perra, M. Karsai, and A. Vespignani, Physical reviewletters , 118702 (2014).[52] M. Keeling and P. Rohani, Modeling Infectious Disease in Hu-mans and Animals (Princeton University Press, 2008).[53] M. Bogu˜na, C. Castellano, and R. Pastor-Satorras, PhysicalReview Letter , 068701 (2013).
SUPPLEMENTARY MATERIAL
Here, we provide supplemental information about the math-ematical derivation and present the sensitivity analysis of theresults to the variation of the main parameters.
ANALYTICAL DERIVATIONS
We consider two types of viruses. The first does not haveaccess to the list of contacts of each victim, thus spread onlyvia the connections activated at each time step t . The sec-ond instead, is able to access the list of contacts in the last τ time-steps. As we will see in details later, at the mean-fieldlevel, the structure of the equations regulating the variation ofthe number of infected individuals at early times is the same.However, for simplicity, let us consider first the first type ofvirus.Each node is assigned to one of Q categories x ∈ [1 , . . . , Q ] that distinguish nodes according the their gullibility λ x , whichdescribes the probability that they will not recognize the threatfor example clicking on the added piece of malicious content,and the time they need to recover for a successful attack, µ − x ( µ x is the recovery rate). Nodes are also characterized by theiractivity a which describes their propensity to engage in socialinteractions in the unit time. To account for observations inreal systems, we extract activity from a heterogenous distribu-tion in particular a power-law F ( a ) = Ba − α with a ∈ [ (cid:15), toavoid divergences. Each active node creates m random con-nections. The target of each communication act is selectedwith probability p within nodes in the same category and with probability − p with nodes in other categories. In both casesthe actual target is selected at random. The virus will be ableto spread only via the connections created by active and in-fected users. In particular, suppose that node i has been com-promised. At time t the node activates and sends m legitimatemessages to m users. During the same time-step, the virus,unbeknownst to i , will send a message with malicious contentto all m users. In these settings, the variation of number ofinfected nodes in each activity class a and category x can bewritten as: d t I xa = − µ x I xa + λ x mS xa p (cid:90) da (cid:48) a (cid:48) I xa (cid:48) N x + (1 − p ) (cid:88) y (cid:54) = x (cid:90) da (cid:48) a (cid:48) I ya (cid:48) N − N y . (3)The first term on the right hand side, describes the recoveryprocess. The second term instead describes susceptible nodesthat are connected by active nodes in the same category thatare infected. These nodes get infected with probability λ x andselected with probability p mN x (where N x is the total num-ber of nodes in the category x ). The third term, accounts forthe same process but in which the susceptible node in activ-ity class a receives a message from active and infected nodesin other categories. Each node is selected with probability (1 − p ) mN − N y by active vertices in class y . At early stages ofthe spreading we can assume that the number of infected to bevery small respect to the susceptible thus we can approximate S xa ∼ N xa . This is equivalent to neglect terms of the order of ( I xa ) . We can also define (cid:82) da (cid:48) a (cid:48) I xa (cid:48) = Θ x , thus summingover all activity classes we get: d t I x = − µ x I x + λ x m p Θ x + (1 − p ) (cid:88) y (cid:54) = x N x N − N y Θ y . (4)In order to characterize the behavior of the number of infectedat such early times, we can write, starting from Eq. 3 the equa-tion for each auxiliary function Θ x . In particular, we can mul-tiply both sides of Eq. 3 for a and integrate over all classes ofactivity. Doing so, we obtain: d t Θ x = − µ x Θ x + λ x m p Θ x (cid:90) da aN xa N x + (1 − p ) (cid:88) y (cid:54) = x N x N − N y Θ y (cid:90) da aN xa N x . (5)where we have multiply and divided the third term for N x . Wecan now define F x ( a ) = N xa N x as the distribution of activitiesin the category x , and thus (cid:82) da aN xa N x = (cid:82) daaF x ( a ) = (cid:104) a (cid:105) x is the average activity in the category. Finally, we let’s define c x,y = N x N − N y which is acts as the mixing probability betweencategorie. In these settings we get: d t Θ x = − µ x Θ x + λ x m (cid:104) a (cid:105) x p Θ x + (1 − p ) (cid:88) y (cid:54) = x c x,y Θ y . (6)Thus we have a system of differential equations made of Q equations. In particular, we have two equations for each x in the form: d t I x = − µ x I x + λ x m p Θ x + (1 − p ) (cid:88) y (cid:54) = x c x,y Θ y = g x .d t Θ x = − µ x Θ x + λ x m (cid:104) a (cid:105) x p Θ x + (1 − p ) (cid:88) y (cid:54) = x c x,y Θ y = h x . (7)The conditions for the spreading can be identified by studyingthe eigenvalues of the Jacobian matrix of such system. TheJacobian can be written as follows: J = ∂g ∂I ∂g ∂I . . . ∂g ∂I Q ∂g ∂ Θ ∂g ∂ Θ . . . ∂g ∂ Θ Q ∂g ∂I ∂g ∂I . . . ∂g ∂I Q ∂g ∂ Θ ∂g ∂ Θ . . . ∂g ∂ Θ Q ... ... . . . ... ... ... . . . ... ∂g Q ∂I ∂g Q ∂I . . . ∂g Q ∂I Q ∂g Q ∂ Θ ∂g Q ∂ Θ . . . ∂g Q ∂ Θ Q ∂h ∂I ∂h ∂I . . . ∂h ∂I Q ∂h ∂ Θ ∂h ∂ Θ . . . ∂h ∂ Θ Q ∂h ∂I ∂h ∂I . . . ∂h ∂I Q ∂h ∂ Θ ∂h ∂ Θ . . . ∂h ∂ Θ Q ... ... . . . ... ... ... . . . ... ∂h Q ∂I ∂h Q ∂I . . . ∂h Q ∂I Q ∂h Q ∂ Θ ∂h Q ∂ Θ . . . ∂h Q ∂ Θ Q (8)Substituting the general terms with the actual partial deriva-tives we get: J = − µ . . . pλ m (1 − p ) λ mc , . . . (1 − p ) λ mc ,Q − µ . . . − p ) λ mc , pλ m . . . (1 − p ) λ mc ,Q ... ... . . . ... ... ... . . . ... . . . − µ Q (1 − p ) λ mc Q, (1 − p ) λ mc Q, . . . pλ Q m . . . − µ + pβ (1 − p ) β c , . . . (1 − p ) β c ,Q . . . − p ) β c , − µ + pβ . . . (1 − p ) β c ,Q ... ... . . . ... ... ... . . . ... . . . − p ) β Q c Q, (1 − p ) β c Q, . . . − µ Q + pβ Q (9)where we defined β x = m (cid:104) a (cid:105) x λ x . It is important to noticethe peculiarities of the Jacobian. The first Q × Q block madeof the partial derivatives of the g x functions in the various I x is a diagonal block that features the recovery rates of eachcategory. The second block on the bottom left side is a Q × Q block of all zeros. Indeed the variables I x do not appear in the h x equations. The adjacent block on the right, features in thediagonal the same function − µ x + pβ x . Due these properties, Q eigenvalues are negative and equal to the negative of eachrecovery rate. The largest eigenvalue instead can be written as Λ max = − (cid:88) x µ x + p (cid:88) x β x + Ξ (10)where Ξ is an algebraic term function of all the β x , µ x and c x,y . We focus on Λ max because the virus will be able tospread if and only if the largest eigenvalue is larger than zero.From this observation we obtain the conditions spreading: R = p (cid:80) x β x + Ξ (cid:80) x µ x > (11)where R is the reproductive number defined as the number ofinfected nodes generated by an initial seed in a fully suscepti-ble population. It is important to mention that for any numberof categories Ξ has an analytical expression. However, sinceit derives from the characteristic equation of the Jacobian ma-trix, Ξ gets more and more complicated as the dimensionality of the matrix increases. Generally speaking for Q categories Ξ is a polynomial of order Q in all variables. Q=1
In case of single category the expression of R becomes: R = βµ (12)In fact, in this limit p = 1 and the Jacobian matrix reduces to J = (cid:20) − µ − µ + β (cid:21) (13)The two eigenvalues are − µ and − µ + β . Thus the diseasewill be able to spread only if β > µ . Q=2
In the case of two categories, Q = 2 , the Jacobian becomes: J = − µ pλ m (1 − p ) λ mc , − µ (1 − p ) λ mc , pλ m − µ + pβ (1 − p ) β c , − p ) β c , − µ + pβ (14)In these settings we have: Ξ = ( µ − µ ) + p ( β − β ) + 2 p ( µ − µ )( β − β ) + 4 β β c , c , ( p − (15)It is important to notice how with two categories, indepen- dently of their sizes c , = c , = 1 . In fact, the two sizes areconstrained by N = N + N . Thus we have: c , = N N − N = N N − N + N = c , = N N − N = N N − N + N = 1 (16)The expression of Ξ reduces to: Ξ = ( µ − µ ) + p ( β + β ) + 2 p ( µ − µ )( β − β ) + 4 β β (1 − p ) (17) Q > As mentioned above, in the most general case of Q cate-gories, the expression of Ξ , becomes quite complex. How-ever, its expression is set unequivocally by the characteristicequation of the Jacobian matrix and can be easily obtainedwith any programming language that allows symbolic com-putations such as Mathematica. The problem can be signifi-cantly simplified in case some of the variables describing thesystem are set. For example in the case of Q = 3 one mightwonder what is the critical value of λ in a system in which β y (with y = [2 , ) and µ y with ( y = [1 , , ) are set. Inthese settings, as shown later on, it is extremely easy to com-pute the largest eigenvalue of the Jacobian for the particularsystem under consideration as function of λ . τ > We now turn the attention to the second type of virus thatis able to access not only the connections establish at time t but also those in previous τ time steps. In order to charac-terize the conditions for the spreading in this case, let us firstunderstand how many people del virus will be able to reachfrom each node of activity a . This number is equal to theout-degree of those nodes. In the case considered in the pre-vious sections τ = 1 , thus the virus was able to reach only thenodes contacted by each active and infected node within thetime-step t . By construction, the out-degree of such nodes is k out ( a ) = ma , since their are active with probability a andwhen active they create m random connections. What aboutfor τ = 2 ?Active nodes at time t might either have m con-nections or m . The first group describes nodes that were notactive at time t − but they were active at time t . The sec-ond group instead describe nodes that were active in both time steps. Thus: k out ( a ) = (1 − a ) am + 2 ma = m ( a + a ) . (18)In fact, nodes of activity a are not active with probability − a and are active two times in a row with probability a (sincethe events are independent). The same reasoning applies for τ = 3 . Here we could have three groups having either degree m , m , and m . As before, the first group describes nodesthat were not active at time t − and t − but they wereactive at time t . The second group instead accounts for all thenodes that were active two times. Finally the third those thatwere active three times. Thus we get: k out ( a ) = ma (1 − a ) + 4 ma (1 − a ) + 3 ma = m ( a + 2 a ) (19)In the case τ = 4 instead we have: k out ( a ) = ma (1 − a ) + 6 ma (1 − a ) + (20) + 9 ma (1 − a ) + 4 ma = m ( a + 3 a ) It is clear that the structure of the out-degree for a general τ can be written as: k out ( a ) = m (cid:2) a + ( τ − a (cid:3) . (21)Within a mean-field approximation, we can approximate theprocess assuming that the virus will try to infected k out ( a ) other nodes as for the case τ = 1 . This is an approximationbecause each active node, at time t , as a quenched list of con-tacts, those established in the time-steps before. The node willnot re-draw them ex novo as in the case τ = 1 . Thus, we canexpect the approximation to be closer to the actual process forsmall values of τ . Within such approach, the structure of theequation is the same as those above, the only different is in the β s since we will have m (cid:104) a (cid:105) x → m (cid:2) (cid:104) a (cid:105) x + ( τ − (cid:104) a (cid:105) x (cid:3) . FEATURES OF THE PHASE SPACE FOR Q = 2 Let’s first consider the case in which µ = µ = µ . Theexpression of R reduces to: R = p ( β + β ) + (cid:112) p ( β + β ) + 4 β β (1 − p )2 µ (22)In the limit p = 0 , nodes in each category will connect justwith nodes in the other. The expression of R becomes: R = √ β β /µ . In the opposite limit, p = 1 , nodes in thetwo categories are separated. Thus we have two independentconditions that have the same mathematical form we encoun-tered for Q = 1 . In fact, we have R = β /µ and R = β /µ .The virus will be able to spread in the system in case either ofthe R x are larger than one. Of course, in case both are largerthan one each group will experience the virus. What happensin case < p < ? It is interesting to notice how the value of R for a general p is bounded by the R x of the two categoriestaken in isolation: min x R x ≤ R ( p ) ≤ max x R x . Before themathematical proof, let us try to develop the intuition behind.Suppose that β > β . Any value of p < , will reduce thespreading power of nodes in the first category. In fact, nodesin category one will be connected to some nodes in category two that are less gullible, or less active, or create a smallernumber of connection (remember that β x = m (cid:104) a (cid:105) x λ x ). Con-versely, nodes in category two, will get in contact with nodesthat increase the spreading potential of the virus. In order toprove this, let us consider the case β > β . We have to showhow R > R ( p ) and R < R ( p ) . Let us consider the firstcondition: β µ > p ( β + β ) + (cid:112) p ( β + β ) + 4 β β (1 − p )2 µ , (23)which is equivalent to: β (2 − p ) − pβ > (cid:112) p ( β + β ) + 4 β β (1 − p ) (24)This condition is respected in case β (2 − p ) − pβ > , p ( β + β ) + 4 β β (1 − p ) > and ( β (2 − p ) − pβ ) >p ( β + β ) + 4 β β (1 − p ) . The first condition implies β > pβ − p , which is always true since β > β was the ini-tial assumption. Furthermore, it is easy to show that equation p ( β + β ) + 4 β β (1 − p ) = 0 as no solution in p , thusthe condition is always respected. Finally, the third conditionimplies β + p β − β p + p β − p (2 − p ) β β > p β + p β + 2 p β β + 4 β β − pβ β (25)that reduces to β > β . The three conditions prove Eq. 23for all p . We have now to prove β µ < p ( β + β ) + (cid:112) p ( β + β ) + 4 β β (1 − p )2 µ , (26)which is equivalent to: β (2 − p ) − pβ < (cid:112) p ( β + β ) + 4 β β (1 − p ) (27)This condition is respected in region in which β (2 − p ) − pβ ≥ , ( β (2 − p ) − pβ ) < p ( β + β ) +4 β β (1 − p ) and p ( β + β ) + 4 β β (1 − p ) ≥ , β (2 − p ) − pβ < .The first two conditions are respected when in the region pβ − p ≤ β < β . The other two instead in the region β < pβ − p . Overall, Eq. 26 is valid in the union of these twothat implies β < β which is exactly the initial assumption.Let’s consider now the general case in which also the tworecovery rates are different. In the limit p = 0 , we have R = √ ( µ − µ ) − β β µ + µ . In the opposite limit instead, p = 1 , the two categories are independent thus we have two condi-tions as before: R = β /µ and R = β /µ . It is inter-esting to notice how in case the two recovery rates are not thesame, the phase space of the process becomes significantlymore complex. In fact, differences in the rate at which nodesrecovers might create interesting non-linear behaviors. In par-ticular, consider a scenario in which the first category featuresa larger β and µ respect to the second. Thus, such nodes aremore prone to infection but recover faster. In case p < , thecoupling between the two categories might boost the spread-ing of the virus, since the node in category one are able toinfect those in two which, although less prone to the diseasestay infected for longer. For a given configuration of param-eters (i.e. setting β s and µ s) we can analytically determinethe value of p above which this phenomenon is observed. Inparticular, let’s assume that β /µ < β /µ . Next, we needto compute the value of p (if any), for which β /µ < R ( p ) .This implies: β µ < p ( β + β ) + Ξ µ + µ (28)that can be written as: β ( µ + µ ) − µ p ( β + β ) < µ (cid:112) ( µ − µ ) + p ( β + β ) + 2 p ( µ + µ )( β − β ) + 4 β β (1 − p ) (29)0It is important to notice how this inequality is at the first orderin p . Indeed, all second order terms cancel out. The value of p that verifies the above inequality lays in the union of two sys-tems of inequalities: i) β ( µ + µ ) − µ p ( β + β ) < andthe quantity inside the square root is larger equal than zero,ii) β ( µ + µ ) − µ p ( β + β ) > , and ( β ( µ + µ ) − µ p ( β + β )) < µ Ξ . Extensive numerical computations show that the values inside the square roots are always posi-tive. Furthermore, the first condition in the first system resultin values of p always larger than one. Thus, the first systemdoes not provide any physical ( p < ) condition. Conversely,the first condition in the second system implies p < whilethe second: p > p ∗ = β ( µ + µ ) − µ ( µ − µ ) − β β µ µ β ( µ + µ )( β + β ) + 2 µ ( µ − µ )( β − β ) − β β µ . (30)Thus, this is the only physical condition necessary to observea reproductive number larger than in each category in isola- tion. Clearly, in the case β /µ < β /µ the condition abovebecomes: p > p ∗ = β ( µ + µ ) − µ ( µ − µ ) − β β µ µ β ( µ + µ )( β + β ) + 2 µ ( µ − µ )( β − β ) − β β µ . (31) p R A) p B) p C) p D) FIG. 4. R as function of p . The shaded area describe the regionin which min x β x /µ x ≤ R ≤ max x β x /µ x . The vertical linedescribe the value of p ∗ from conditions Eq. 30 and Eq. 31. In panelsA-B we set µ = 10 − , µ = 5 × − , m = 4 , λ = 0 . , λ = 0 . (A) and λ = 0 . (B). In panels C-D we set µ = 5 × − , µ =3 × − , m = 4 , λ = 0 . , λ = 0 . (C) and λ = 0 . (D). In Figure 4 we verify the above condition. In particular, we setthe values of β x and µ x and plot R from Eq. 11 as functionof p . In particular, we consider that nodes are assigned tothe categories randomly. The shaded area is the region where min x β x /µ x ≤ R ≤ max x β x /µ x . The vertical line showthe value of p ∗ determined from the condition derived above.It is clear how for a given setting, there might be a value of p above which the reproductive number gets indeed larger thanthe the R x of each category in isolation.It is important to stress how the region of the phase spacein which we observe this phenomenon is generally speakingquite limited. In fact, it might happen only in case the categorywith the larger recovery rates has also the larger gullibility. InFigure 5 we show as contour plots the region of parameterswhere the reproductive number of the system is larger thanthat correspondent value in the two categories in isolation. Inparticular, we set λ = 0 . , λ = 0 . (A), λ = 0 . (B), λ = 0 . (C), λ = 0 . (D) and show as function of µ and A) B) C) D) FIG. 5. We show as function of µ and µ the region of parametersin which the reproductive number of system is larger than the corre-spondent values computed in each category in isolation. The colorsrefer to the value of p (calculated from Eq. 30 and Eq. 31) abovewhich this phenomenon is observed. We set λ = 0 . , λ = 0 . (A), λ = 0 . (B), λ = 0 . (C), λ = 0 . (D) µ the value of p ∗ . It is clear this region increases as the dif-ference between the two gullibilities increases. It is importantto notice how the expression for p ∗ is perfectly in line with thecase in which µ = µ . Indeed, in this limit we get p ∗ > which implies, as expected, that the necessary condition tohave a reproductive number larger than in each category inisolation is to have different recovery rates. NUMERICAL SIMULATIONS
In this section we will present the sensitivity analysis tothe model’s parameters. We will first consider two categories( Q = 2 ). As shown in the main text, we adopted two main ap-proaches to assign node to categories. The first is at random,the second is instead in decreasing order to activity. In partic-ular, we order activity in decreasing order and then assign thefirst gN nodes to the first category and the remaining to the1 L A) B) C) D) E) F) FIG. 6. Lifetime of the SIS process (A-C) and contour plot of R (D-F). In A-B-D-E nodes are randomly assigned to two categories,in C-F instead in decreasing order of activity. We set p = 0 . (A-D), p = 0 . (B-C-E-F). In A-C we fix N = 2 × , m = 4 , α = − . , µ = µ = 10 − , λ = 0 . , Y = 0 . , and . of randominitial seeds. We plot the median and confidence intervals in simulations per point. The solid lines come from Eq. 11, and thedashed lines are the analytical threshold in case of a single category.In the contour plot we set µ = µ = 10 − . second. Thus (cid:104) a (cid:105) = (cid:82) a c daaF ( a ) and (cid:104) a (cid:105) = (cid:82) a c (cid:15) daaF ( a ) and a c is determined in such a way that the fraction of nodesin the first class is g . This can be easily done imposing: (cid:90) a c F ( a ) da = g. (32)Since F ( a ) = − α − (cid:15) − α a − α we get: a c = (cid:2) − g (1 − (cid:15) − α ) (cid:3) − α (33)It is important to notice that in Eq. 11 the expression of (cid:104) a (cid:105) x in the two assignment scenarios is slightly different. In par-ticular we defined (cid:104) a (cid:105) x = (cid:82) daF x ( a ) a = (cid:82) da N xa N x a . In casenodes are assigned randomly to the two categories we havethat F x ( a ) ∼ F ( a ) since N xa = N a /g and N x = N/g (where g is the fraction of node in the general category x in this case).Thus, (cid:104) a (cid:105) x = (cid:104) a (cid:105) for the two categories. In case instead nodesare assigned in decreasing order of activity (cid:104) a (cid:105) x = (cid:104) a (cid:105) /g . Infact, in this limit N xa = N a (since nodes are assigned to cate-gories as function of their activity) but N x = gN .In Figure 6, we consider the case of µ = µ in case of ran-domly selected nodes for p = 0 . (A-D) and p = 0 . (B-E).We also consider the case of correlation between category andactivity for p = 0 . (C-F). Respect to Figure 1 of the maintext, we used a different value of λ = 0 . . Across the boardthe analytical solution match perfectly the numerical estima-tion of the threshold.In Figure 7 we consider the general case of different recoveryrates. In particular we consider a different set of values respectto those in the main text. In particular, we set µ = 10 − , µ = 10 − , λ = 0 . and m = 4 . In panels A-B-D-E weconsider random assignment of nodes to categories. In C-Fwe consider the correlation between activity and category. Weassign to category one to most active nodes. Also, in panelsB-C-E-F we considered p = 0 . while in panel A-D we set p = 0 . . Overall, the figure confirms the validity of the the-oretical approach and highlights one more time the effects of c L A) c B) c C) D) E) F) FIG. 7. Lifetime of the SIS process (A-C) and contour plot of R (D-F). In A-B-D-E nodes are randomly assigned to two categories,in C-F instead in decreasing order of activity. We set p = 0 . (A-D), p = 0 . (B-C-E-F). In A-C we fix N = 2 × , m = 4 , α = − . , µ = 10 − , µ = 10 − , λ = 0 . , Y = 0 . , and . of randominitial seeds. We plot the median and confidence intervals in simulations per point. The solid lines come from Eq. 11, and thedashed lines are the analytical threshold in case of a single category.In the contour plot we set λ = 0 . and λ = 0 . . c L A) c B) c C) D) E) F) FIG. 8. Lifetime of the SIS process (A-C) and contour plot of R (D-F). In A-B-D-E nodes are randomly assigned to two categories,in C-F instead in decreasing order of activity. We set p = 0 . (A-D), p = 0 . (B-C-E-F). In A-C we fix N = 2 × , m = 6 , α = − . , µ = 10 − , µ = 10 − , λ = 0 . , Y = 0 . , and . of randominitial seeds. We plot the median and confidence intervals in simulations per point. The solid lines come from Eq. 11, and thedashed lines are the analytical threshold in case of a single category.In the contour plot we set λ = 0 . and λ = 0 . . correlations between category assignment and activity that re-duce the non-active phase space (see panel F). Furthermore,it is important to notice how the critical value in case of asingle category with a recovery rate average of the two herewould be λ c = 2 . (not shown in the figure) which impliesthat the virus would not be able to spread since all the gulli-bilities should be smaller or equal to . This confirms the im-portance of accounting for the presence of different categoriesof users in order to correctly capture the spreading power ofthe virus.In Figure 8 we test the sensitivity to the parameter m . In themain text as well as in many of the other plots we set m = 4 .Here, we fix instead m = 6 keeping all the other parame-ters the same as in the Figure 7. The analytical solutions onemore time match the numerical simulations and the contourplots confirm the picture discussed in the main text and allthe other similar plots. In Figure 9 we test the sensitivity tothe exponent of the activity distribution. In all the other plotswe set α = − . , here instead we consider α = − . . Weconsidered a scenario in which the recovery rates of the two2 c L A) c B) c C) D) E) F) FIG. 9. Lifetime of the SIS process (A-C) and contour plot of R (D-F). In A-B-D-E nodes are randomly assigned to two categories,in C-F instead in decreasing order of activity. We set p = 0 . (A-D), p = 0 . (B-C-E-F). In A-C we fix N = 2 × , m = 6 , α = − . , µ = 10 − , µ = 10 − , λ = 0 . , Y = 0 . , and . of randominitial seeds. We plot the median and confidence intervals in simulations per point. The solid lines come from Eq. 11, and thedashed lines are the analytical threshold in case of a single category.In the contour plot we set λ = 0 . and λ = 0 . . categories is the same, set λ = 0 . , m = 6 and consider two different values of p . As clear from the figure, also in thiscase the analytical estimation matches the numerical simula-tions. Furthermore, it is interesting to notice how, in case offaster decay of the activity distribution (i.e. smaller value ofthe exponent α ), the threshold of the correlated case (panel C-F) is closer to the scenario of a single category (dashed line).Indeed, the average activity of the more active category getscloser to the average activity of the whole network. Q = 3 Here we consider the case of three categories. For sim-plicity let’s consider nodes are assigned to the categories atrandom and that categories have the same size N x = N/ .Also, let’s us set the values of β x with x = [2 , , µ x with x = [1 , , , m = 4 , and assume that links are created ran-domly between categories thus p = 1 / . In particular, if weset β = β = 0 . , µ = µ = µ = 0 . , we can use Eq. 9to obtain the critical value of λ . In particular, the generalexpression of the Jacobian is: J = − µ pλ m (1 − p ) λ mc , (1 − p ) λ mc , − µ − p ) λ mc , pλ m (1 − p ) λ mc , − µ (1 − p ) λ mc , (1 − p ) λ mc , pλ m − µ + pβ (1 − p ) β c , (1 − p ) β c , − p ) β c , − µ + pβ (1 − p ) β c , − p ) β c , (1 − p ) β c , − µ + pβ (34)Since the categories have the same size: c x,y = N x N − N y = N N − N = 12 (35)Plugging all the values and solving for λ we obtain: λ c = 4255 (36) In Figure 10 we show the comparison between the analyti-cal prediction and the numerical simulations which perfectlymatches.3 c L FIG. 10. We show the lifetime of the SIS process in case of Q = 3 as function of λ . The vertical line describes the analytical estima-tion of its critical value. In the simulation we set β = β = 0 . , µ = µ = µ = 0 . , N = 3 × , m = 4 , α = − . , (cid:15) = 10 − and run simulations for each data point. We show the50%