The spreading of computer viruses on time-varying networks

TL;DR

This paper presents a theoretical framework for understanding how computer viruses spread on time-varying social networks, considering user susceptibility and network dynamics, revealing increased system fragility due to temporal interactions.

Contribution

It introduces a novel model that incorporates temporal network properties and user heterogeneity, providing analytical insights into virus invasion thresholds.

Findings

Temporal coupling of user categories increases system fragility.

Network dynamics significantly influence virus spreading.

Analytical invasion thresholds derived for realistic scenarios.

Abstract

Social networks are the prime channel for the spreading of computer viruses. Yet the study of their propagation 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 viruses propagating on temporal networks featuring Q categories of susceptibility and derive analytically the invasion threshold. We found that the temporal coupling of categories might increase the fragility of the system to cyber threats. Our results show that networks' dynamics and their interplay with users features are crucial for the spreading of computer viruses.

Figures

Equations

$$d_t{I}_a^x=-\mu I_a^x+{\lambda }_{x}m{S}_a^x\times$$

$$\left[p\int d{a}^{\prime }{a}^{\prime }\frac{I_{a^{\prime }}^x}{N^x}+(1-p)\sum _{y\ne x}\int d{a}^{\prime }{a}^{\prime }\frac{I_{a^{\prime }}^y}{N-N^y}\right].$$

$$d_t{I}^x=-{\mu }_x{I}^x+{\lambda }_{x}m\left[p{\theta }^x+(1-p)N^x\sum _{y\ne x}{\theta }^y/(N-N^y)\right],$$

$$d_t{\theta }^x=-{\mu }_x{\theta }^x+$$

$$m{\lambda }_x⟨a{⟩}_x\left[p{\theta }^x+(1-p)N^x\sum _{y\ne x}{\theta }^y/(N-N^y)\right].$$

$$R_0=\frac{p{\sum }_x{\beta }_x+\Xi }{{\sum }_x{\mu }_x}>1,$$

$$d_t{I}_a^x=-{\mu }_x{I}_a^x+{\lambda }_{x}m{S}_a^x\left[p\int d{a}^{\prime }{a}^{\prime }\frac{I_{a^{\prime }}^x}{N^x}+(1-p)\sum _{y\ne x}\int d{a}^{\prime }{a}^{\prime }\frac{I_{a^{\prime }}^y}{N-N^y}\right].$$

$$d_t{I}^x=-{\mu }_x{I}^x+{\lambda }_{x}m\left[p{\Theta }^x+(1-p)\sum _{y\ne x}\frac{N^x}{N-N^y}{\Theta }^y\right].$$

$$d_t{\Theta }^x=-{\mu }_x{\Theta }^x+{\lambda }_{x}m\left[p{\Theta }^x\int da\frac{a{N}_a^x}{N^x}+(1-p)\sum _{y\ne x}\frac{N^x}{N-N^y}{\Theta }^y\int da\frac{a{N}_a^x}{N^x}\right].$$

$$d_t{\Theta }^x=-{\mu }_x{\Theta }^x+{\lambda }_{x}m⟨a{⟩}_x\left[p{\Theta }^x+(1-p)\sum _{y\ne x}{c}_{x,y}{\Theta }^y\right].$$

$$d_t{I}^x$$

$$d_t{\Theta }^x$$

$$J=\left[\begin{array}{cccccccc}\frac{\partial g^1}{\partial I^1}& \frac{\partial g^1}{\partial I^2}& \dots & \frac{\partial g^1}{\partial I^Q}& \frac{\partial g^1}{\partial {\Theta }^1}& \frac{\partial g^1}{\partial {\Theta }^2}& \dots & \frac{\partial g^1}{\partial {\Theta }^Q}\\ \frac{\partial g^2}{\partial I^1}& \frac{\partial g^2}{\partial I^2}& \dots & \frac{\partial g^2}{\partial I^Q}& \frac{\partial g^2}{\partial {\Theta }^1}& \frac{\partial g^2}{\partial {\Theta }^2}& \dots & \frac{\partial g^2}{\partial {\Theta }^Q}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ \frac{\partial g^Q}{\partial I^1}& \frac{\partial g^Q}{\partial I^2}& \dots & \frac{\partial g^Q}{\partial I^Q}& \frac{\partial g^Q}{\partial {\Theta }^1}& \frac{\partial g^Q}{\partial {\Theta }^2}& \dots & \frac{\partial g^Q}{\partial {\Theta }^Q}\\ \frac{\partial h^1}{\partial I^1}& \frac{\partial h^1}{\partial I^2}& \dots & \frac{\partial h^1}{\partial I^Q}& \frac{\partial h^1}{\partial {\Theta }^1}& \frac{\partial h^1}{\partial {\Theta }^2}& \dots & \frac{\partial h^1}{\partial {\Theta }^Q}\\ \frac{\partial h^2}{\partial I^1}& \frac{\partial h^2}{\partial I^2}& \dots & \frac{\partial h^2}{\partial I^Q}& \frac{\partial h^2}{\partial {\Theta }^1}& \frac{\partial h^2}{\partial {\Theta }^2}& \dots & \frac{\partial h^2}{\partial {\Theta }^Q}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ \frac{\partial h^Q}{\partial I^1}& \frac{\partial h^Q}{\partial I^2}& \dots & \frac{\partial h^Q}{\partial I^Q}& \frac{\partial h^Q}{\partial {\Theta }^1}& \frac{\partial h^Q}{\partial {\Theta }^2}& \dots & \frac{\partial h^Q}{\partial {\Theta }^Q}\end{array}\right]$$

$$J=\left[\begin{array}{cccccccc}-{\mu }_1& 0& \dots & 0& p{\lambda }_{1}m& (1-p){\lambda }_{1}m{c}_{1,2}& \dots & (1-p){\lambda }_{1}m{c}_{1,Q}\\ 0& -{\mu }_2& \dots & 0& (1-p){\lambda }_{2}m{c}_{2,1}& p{\lambda }_{2}m& \dots & (1-p){\lambda }_{2}m{c}_{2,Q}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & -{\mu }_Q& (1-p){\lambda }_{2}m{c}_{Q,1}& (1-p){\lambda }_{2}m{c}_{Q,2}& \dots & p{\lambda }_{Q}m\\ 0& 0& \dots & 0& -{\mu }_1+p{\beta }_1& (1-p){\beta }_1{c}_{1,2}& \dots & (1-p){\beta }_1{c}_{1,Q}\\ 0& 0& \dots & 0& (1-p){\beta }_2{c}_{2,1}& -{\mu }_2+p{\beta }_2& \dots & (1-p){\beta }_2{c}_{2,Q}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 0& (1-p){\beta }_Q{c}_{Q,1}& (1-p){\beta }_2{c}_{Q,2}& \dots & -{\mu }_Q+p{\beta }_Q\end{array}\right]$$

$${\Lambda }_{max}=-\sum _x{\mu }_x+p\sum _x{\beta }_x+\Xi$$

$$R_0=\frac{p{\sum }_x{\beta }_x+\Xi }{{\sum }_x{\mu }_x}>1$$

$$R_0=\frac{\beta }{\mu }$$

$$J=\left[\begin{array}{cc}-\mu & 0\\ 0& -\mu +\beta \end{array}\right]$$

$$J=\left[\begin{array}{cccc}-{\mu }_1& 0& p{\lambda }_{1}m& (1-p){\lambda }_{1}m{c}_{1,2}\\ 0& -{\mu }_2& (1-p){\lambda }_{2}m{c}_{2,1}& p{\lambda }_{2}m\\ 0& 0& -{\mu }_1+p{\beta }_1& (1-p){\beta }_1{c}_{1,2}\\ 0& 0& (1-p){\beta }_2{c}_{2,1}& -{\mu }_2+p{\beta }_2\end{array}\right]$$

$${\Xi }^2=({\mu }_1-{\mu }_2{)}^2+p^2({\beta }_1-{\beta }_2{)}^2+2p({\mu }_2-{\mu }_1)({\beta }_1-{\beta }_2)+4{\beta }_1{\beta }_2{c}_{1,2}{c}_{2,1}(p-1{)}^2$$

$$c_{1,2}=\frac{N^1}{N-N^2}=\frac{N^1}{N-N+N^1}=c_{2,1}=\frac{N^2}{N-N^1}=\frac{N^2}{N-N+N^2}=1$$

$${\Xi }^2=({\mu }_1-{\mu }_2{)}^2+p^2({\beta }_1+{\beta }_2{)}^2+2p({\mu }_2-{\mu }_1)({\beta }_1-{\beta }_2)+4{\beta }_1{\beta }_2(1-2p)$$

$$k^{out}(a)=(1-a)am+2m{a}^2=m(a+a^2).$$

$$k^{out}(a)=ma(1-a{)}^2+4m{a}^2(1-a)+3m{a}^3=m(a+2a^2)$$

$$k^{out}(a)$$

$$k^{out}(a)=m\left[a+(\tau -1)a^2\right].$$

$$R_0=\frac{p({\beta }_1+{\beta }_2)+\sqrt{p^2({\beta }_1+{\beta }_2{)}^2+4{\beta }_1{\beta }_2(1-2p)}}{2\mu }$$

$$\frac{{\beta }_1}{\mu }>\frac{p({\beta }_1+{\beta }_2)+\sqrt{p^2({\beta }_1+{\beta }_2{)}^2+4{\beta }_1{\beta }_2(1-2p)}}{2\mu },$$

$${\beta }_1(2-p)-p{\beta }_2>\sqrt{p^2({\beta }_1+{\beta }_2{)}^2+4{\beta }_1{\beta }_2(1-2p)}$$

$$4{\beta }_1^2+p^2{\beta }_1^2-4{\beta }_1^{2}p+p^2{\beta }_2^2-2p(2-p){\beta }_1{\beta }_2>p^2{\beta }_1^2+p^2{\beta }_2+2p^2{\beta }_1{\beta }_2+4{\beta }_1{\beta }_2-8p{\beta }_1{\beta }_2$$

$$\frac{{\beta }_2}{\mu }<\frac{p({\beta }_1+{\beta }_2)+\sqrt{p^2({\beta }_1+{\beta }_2{)}^2+4{\beta }_1{\beta }_2(1-2p)}}{2\mu },$$

$${\beta }_2(2-p)-p{\beta }_1<\sqrt{p^2({\beta }_1+{\beta }_2{)}^2+4{\beta }_1{\beta }_2(1-2p)}$$

Full text

Alongside clear societal and economic benefits, modern technology exposes us to serious challenges. In particular, the spreading of malicious content online, often based on ingenious deception strategies, is one of the most pressing because it poses serious threats to our privacy, finances, and safety Kayes and Iamnitchi (2017). Victims of a typical social engineering attack Heartfield and Loukas (2016) may receive a message containing a malicious link or file, appearing to originate from a friend or other trusted entity. If opened, it may compromise the computer, access personal information, and spread the virus further unbeknownst to the victim. Recent research has shown how the susceptibility of individuals to such attacks is not homogenous and depends on several features such as age, prior training, computer proficiency, familiarity with social network platforms, among others Heartfield and Loukas (2018); Heartfield et al. (2017, 2016). Furthermore, the properties of real networks are known to facilitate the propagation of such processes Lloyd and May (2001); Balthrop et al. (2004); Pastor-Satorras and Vespignani (2001); Moreno and Vázquez (2003); Newman (2002, 2010); Pastor-Satorras et al. (2015); Barrat et al. (2008); Yang et al. (2013); Yang and Yang (2014). In particular, the heterogeneity in contact patterns makes socio-technical systems quite fragile to biological and digital threats.

The study of these phenomena has largely neglected the complex temporal nature of online contact patterns in favor of static and time-aggregated approaches Holme (2015); Holme and Saramäki (2012). These approximations might be fitting. Indeed, in the past, computer viruses would spread mainly via email networks, targeting the address books of victims, which contain contacts lists Newman et al. (2002). However, not many people create such lists any more and access to them is restricted Balthrop et al. (2004). In the context of social or biological contagions, neglecting the temporal nature of the networks where the processes unfold has been shown to induce misrepresentations of their spreading potential. In fact, the order and concurrency of connections is key Barrat and Cattuto (2015); Perra et al. (2012a, b); Ribeiro et al. (2013); Liu et al. (2014a, 2013); Ren and Wang (2014); Starnini et al. (2013, 2012); Valdano et al. (2015); Scholtes et al. (2014); Williams and Musolesi (2016); Rocha and Masuda (2014); Takaguchi et al. (2012a); Rocha and Blondel (2013); Ghoshal and Holme (2006); Sun et al. (2015); Mistry et al. (2015); Pfitzner et al. (2013); Takaguchi et al. (2012b, 2013); Holme and Liljeros (2014); Holme and Masuda (2015); Wang et al. (2016); Gonçalves and Perra (2015). To the best of our knowledge, beside some early work on the spreading of viruses via Bluetooth among mobile phones Wang et al. (2009), the study of the propagation of cyber threats considering the temporal nature of social interactions is still missing. Furthermore, with few exceptions Peng et al. (2017), the literature devoted to the study of computer viruses unfolding on networks typically neglects that the susceptibility of online users is not homogenous. Conversely, the literature that studies the susceptibility of users to cyber threats traditionally focuses on single users neglecting their connections.

To tackle these limitations, here we introduce a theoretical framework to study the spreading of computer viruses, based on social engineering deception strategies, on time-varying networks. We model users’ interactions using a time-varying network model and consider two types of viruses. The first mimics threats that can propagate only via connections activated at each time step. The second, on the contrary, considers viruses able to access also information about past connections. We investigate the impact of different classes of susceptibility considering that they might also influence the link formation process. In all cases, we analytically derive the conditions regulating the spreading of the virus. Interestingly, these are defined by the interplay between the features of the cyber threats, the categories of susceptibility and their time-varying connectivity. Furthermore, in some scenarios, the coupling between categories creates a complex phenomenology that favors the spreading of the virus. These results have the potential to initiate future efforts aimed at describing more realistically the spreading of computer viruses on online social networks.

We consider a population of $N$ online users which exchange messages in a time-varying network. Nodes are assigned to one of $Q$ categories describing their susceptibility to cyber threats measured in terms of their gullibility and time needed to recover from successful attacks. Since susceptibility is linked to demographic features, we consider that the membership to a category might influence the link creation process. In fact, homophily is a strong social mechanism known to affect the structure and organization of ties McPherson et al. (2001). We model the contact patterns between users with a generalization of the activity-driven framework Perra et al. (2012b); Karsai et al. (2014); Ubaldi et al. (2016); Tizzani et al. (2018). Here, nodes feature an activity $a$ describing their propensity to initiate communications. Activities are extracted from a distribution $F(a)$ which, as observations in real systems have shown, is typically heterogenous Perra et al. (2012b); Ribeiro et al. (2013); Ubaldi et al. (2016); Tomasello et al. (2014). We select power-law distributions $F(a)\sim a^{-\alpha}$ with $a\in[\epsilon,1]$ to avoid divergences. At each time step nodes are active with probability $a\Delta t$. Active nodes select $m$ others and create directed (out-going) links which mimic messages.

In the simplest version of activity-driven networks the selection is random and memoryless Perra et al. (2012b). Here, we propose a variation: with probability $p$ each target is selected, at random, among the group of nodes in the same category, and with probability $1-p$ among the nodes in any other category. In other words, $p$ tunes the homophily level in the network with respect to susceptibility to cyber threats. At time $t+\Delta t$ all edges are deleted and the process starts from the beginning. Unless specified otherwise, links have a duration $\Delta t$. Without loss of generality we set $\Delta t=1$. The model is clearly a simplification of real interactions. However, it offers simple, yet non trivial, settings to study the effects of temporal connectivity patterns on contagion processes unfolding at a comparable time-scale with respect to the evolution of connections Perra et al. (2012b, a); Karsai et al. (2014); Liu et al. (2014b).

We describe the propagation of a computer virus adopting the prototypical SIS model Keeling and Rohani (2008); Barrat et al. (2008). At each time step $t$ the virus, unbeknownst to the victims, sends a message, with malicious content, to all the nodes genuinely contacted at $t$ (virus type 1) or within $t-\tau$ time-step (virus type 2). The focus is not defining the optimal set of nodes to maximize/minimize the damage. Thus, we select randomly a small percentage ($0.5\%$) of nodes as initial seeds. In these settings, susceptible nodes of class $x\in[1,\ldots,Q]$, that receive a malicious message, become infectious with probability $\lambda_{x}$ which defines their gullibility. They recover and become susceptible again with rate $\mu_{x}$. Assuming that nodes with the same value of activity and in the same category are statistically equivalent, we group nodes according to the two features. At each time step, we call $S_{a}^{x}$ and $I_{a}^{x}$ the number of nodes susceptible and infected in activity class $a$ and category $x$. Clearly $\int daS_{a}^{x}=S^{x}$, $\int daI_{a}^{x}=I^{x}$, $\sum_{x}S^{x}=S$, and $\sum_{x}I^{x}=I$. Furthermore, $N_{a}^{x}$ describes the number of nodes of activity $a$ in category $x$, thus $\int daN_{a}^{x}=N^{x}$ and $\sum_{x}N^{x}=N$. In these settings, we can represent the variation of the number of infected nodes of activity $a$ in category $x$ as:

[TABLE]

The first term on the right hand side accounts for the recovery process. The second and third terms capture susceptible nodes that receive messages from active and infected vertices in the same (second) or different (third) category, and get infected as a result. With respect to the typical biological contagion process, here transmission is asymmetric. Only nodes receiving a message from an infected person might be exposed to the virus. Thus, not only the order of connections, but also their direction is a crucial ingredient for the spreading. Since the links are created randomly, each node is selected with a probability $pm/N^{x}$ by nodes in the same category or $(1-p)m/(N-N^{y})$ by nodes in other categories. The total number of nodes is constant thus $S_{a}^{x}=N_{a}^{x}-I_{a}^{x}$ and at the early stages of the spreading we can assume that the number of infected nodes is very small: $S_{a}^{x}\sim N_{a}^{x}$. By integrating across all activities Eq. 3 we get:

[TABLE]

where we define $\theta^{x}=\int daaI_{a}^{x}$. By multiplying both sides of Eq. 3 for $a$ and integrating across all the activities we obtain

[TABLE]

The virus is able to spread, if and only if the largest eigenvalue of the Jacobian matrix of the system of differential equations in $I^{x}$ and $\theta^{x}$ is larger than zero Perra et al. (2012b). As shown in details in the Supplementary Material (SM) this implies:

[TABLE]

where $R_{0}$ is the basic reproductive number defined as the average number of infected nodes generated, in a fully susceptible population, by an infected individual Keeling and Rohani (2008), $\beta_{x}=m\lambda_{x}\langle a\rangle_{x}$ and $\Xi$ is a function of the interplay between the average activation, infection and recovery rate of each category as well as of the mixing between categories.

To understand the dynamics, let us consider a particular case in which the system is characterized by only two categories. Furthermore, let us consider, as first scenario, that all nodes have the same recovery rate. In these settings we have $\Xi^{2}=p^{2}(\beta_{1}+\beta_{2})^{2}+4\beta_{1}\beta_{2}(1-2p)$. The condition for the spreading, even with only two classes, is a non linear function of the average activity of each category, the infection probabilities per contact and the homophily. In the limit $p=0$, nodes in a category connects only with vertices in the other and the expression reduces to $R_{0}=\frac{\sqrt{\beta_{1}\beta_{2}}}{\mu}$. In the limit $p=1$ instead, interactions are only between nodes in the same category. The system is effectively split in two disconnected networks and there are two independent conditions $R_{0}^{x}=\beta_{x}/\mu$. For a general $p$, these two values confine $R_{0}$: $\min_{x}R_{0}^{x}\leq R_{0}(p)\leq\max_{x}R_{0}^{x}$. In fact, any value of $p<1$ will reduce the spreading power of the category characterized with the largest $R_{0}^{x}$ as some connections will be established with nodes where the virus finds it harder to spread (see SM for the proof).

In Fig. 1-A-C, we compare analytical predictions with numerical simulations. We set $\lambda_{2}=0.2$ and use Eq. 2 to estimate the critical value of $\lambda_{1}$ for which $R_{0}\equiv 1$. On the $y$-axis we plot the lifetime of the process defined as the time that the virus needs either to die out or to reach a fraction $Y$ of the population Boguña et al. (2013). The lifetime acts as the susceptibility of a second order phase transition and allows a precise numerical estimation of the threshold of SIS processes Boguña et al. (2013). In panels A-B we consider a scenario in which nodes are assigned randomly to one of the two categories. Thus the average activity in the two is the same and set $p=0.9$ and $p=0.4$ respectively. The analytical value of the threshold (vertical solid line) perfectly matches the numerical estimation. For $p=0.9$ the threshold is smaller than for $p=0.4$ and closer to the threshold of a system with a single category (dashed lines). For smaller values of homophily, instead, the critical conditions are driven by the interplay between the activation rates and gullibility of the two categories. Panels D-E show the analytical value of $R_{0}$ as a function of $\lambda_{1}$ and $\lambda_{2}$ for the two values of $p$. The grey regions 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 will spread if above the threshold in at least one category independently of the other. In the opposite limit, on the contrary, the two categories get intertwined and a small value of the infection probability in one category should be associated to a progressively large value in the other.

In panels C-F we consider that the first category contains a fraction $g$ of nodes selected in decreasing order of activity. Thus, this category contains the $gN$ most active nodes, while the other the $(1-g)N$ least active (see SM). To compare with panel B, we set $g=0.5$ and $p=0.4$. First, the analytical threshold nicely matches the numerical simulations. Second, although the other parameters are the same used in panel B, the critical value of the gullibility of the first class is smaller. Thus, correlations between activity and gullibility facilitate the spreading. This is confirmed in panel F where the active phase space features a region in which the spreading is completely dominated by the category of most active nodes. Overall, all the plots show the importance of distinguishing nodes according to their gullibility. Indeed, neglecting the presence of different classes of users might induce a strong misrepresentation of the virus propagation (dashed lines).

Let us next consider a second scenario where categories differentiate also for the time needed to recover from a successful attack. For two categories, we can write $\Xi^{2}=(\mu_{1}-\mu_{2})^{2}+p^{2}(\beta_{1}+\beta_{2})^{2}+2p(\mu_{2}-\mu_{1})(\beta_{1}-\beta_{2})+4\beta_{1}\beta_{2}(1-2p)$. Interestingly, we have the same terms that appeared in the first scenario, plus two that feature the difference between the recovery rates and $\beta$s of the two categories. Thus $R_{0}$ is a function of the interplay between the activities, gullibilities and recovery rates. In the limit $p=0$, each category only connects with nodes in the other, the two groups are coupled and the threshold reads $R_{0}=\frac{\sqrt{(\mu_{1}-\mu_{2})^{2}+4\beta_{1}\beta_{2}}}{\mu_{1}+\mu_{2}}$. In the limit $p=1$ instead, the two categories are completely de-coupled and the threshold becomes, as before, $R_{0}=\beta_{x}/\mu_{x}$. As shown in Fig. 2-D-H, for a general value of $p$ the reproductive number is not bounded, as before, by the values of $R_{0}^{x}$ computed in the two classes separately (see SM). In Fig. 2-D, we assign nodes randomly to each category, fix $\beta_{x}$ and $\mu_{x}$ and compute $R_{0}$ as a function of $p$. In the shaded area $\min_{x}R_{0}^{x}\leq R_{0}(p)\leq\max_{x}R_{0}^{x}$. Interestingly, after a $p^{*}$ (vertical dashed line), which as shown in the SM can be computed analytically, we enter in a regime where $R_{0}(p)>\max_{x}R_{0}^{x}$. Thus, only specific values of the coupling between categories might induce the virus to spread faster in the combined system than in each single category in isolation. However, this non linear effect is found only in a small fraction of the phase space see Fig. 2-H. The necessary, but not sufficient condition, is that two categories differentiate both for gullibility and recovery rates in such a way that one is more gullible and recovers faster than the other. In this regime, the right mixing between the two might create a feedback loop that makes the system more fragile.

Fig. 2-A-C shows a good match between the analytical (solid vertical lines) and numerical thresholds in case of nodes are assigned at random (A-B) or in decreasing order of activity (C) to the two categories. We fix two different recovery rates, $\lambda_{2}$, and use $\lambda_{1}$ as order parameter. Panels A-B differ in the value of the homophily $p$. We set $p=0.9$ in A, while $p=0.4$ in B-C. The presence of a category of nodes characterized by a smaller value of recovery rate pushes the threshold to smaller values with respect to the first scenarios (Fig. 1). As before, the value of the threshold estimated considering only a single category, characterized by the average recovery rate of the two, (dashed lines) leads to a misrepresentation of the spreading power of the virus, especially for smaller values of homophily (see panel B).

The effect of $p$ on the critical value of $\lambda_{1}$ is similar to the first scenario. In fact, even when categories differentiate by the recovery rates, high values of homophily push the critical point to smaller values. However, here the difference between the two is less significant than in Fig. 1. In Fig. 2-E-F, we show the analytical value of $R_{0}$ as function of $\mu_{1}$ and $\mu_{2}$. Interestingly, the sub-critical region, for $p=0.4$, is smaller than for $p=0.9$. This is in contrast to what was observed in the corresponding plots for the first scenario and highlights once again the complex phenomenology introduced by the interplay of different recovery rates. In Fig. 2-C-G we investigate the effect of correlations. In case that the most active nodes are able to recover quickly from the attack, the virus is able to spread only if the gullibility of such users is higher than in the corresponding case without correlations (panel B). This is confirmed in panel G, where we see that correlations between recovery rates significantly change the active region.

Finally, we turn our attention to a second type of virus able to access also past contacts of infected users within a time window $\tau$. As before, the virus propagates via active infected nodes, but at each time $t$ active users might infect their contacts in a time-window $(t-\tau,t]$. Within a mean-field approximation, we can adopt the same equations described above and change the probability that a node in each activity class receives a message by active and infected nodes. In this case, the out-degree of each active node is not $m$, but a function of $\tau$: $k^{out}(a)=m\left[a+(\tau-1)a^{2}\right]$ (see SM). To grasp the derivation, consider the simplest scenario in which $\tau=2$. In this case, active nodes might have either $m$ or $2m$ contacts in two time steps. The first class describes nodes that are active at time $t$ but were not active at time $t-1$; whereas the second, nodes that were active in both time steps. Thus the out-degree of these nodes, on average, is $k^{out}(a)=ma(1-a)+2ma^{2}$. As shown in the SM, the condition for the spreading has the same structure of Eq. 2 where, however, the value of $\beta$s are changed with the following transformation $m\rightarrow m\left[\langle a\rangle+(\tau-1)\langle a^{2}\rangle\right]$. Thus, the larger the visibility of past connections, from the virus point of view, the larger $R_{0}$. Intuitively this is due to the fact that the virus, for large values of $\tau$, is able to access more contacts, which results in a larger spreading potential. This observation nicely shows how neglecting the temporal nature of connectivity patterns in favor of static (or time integrated) approximations might lead to a poor description of the propagation of viruses that do not have access to contacts lists or past connections. In Fig. 3 we show the comparison between analytical (solid lines) and numerical values of the threshold for different values of $\tau$. To isolate the effect of $\tau$ we considered two categories, a single recovery rate, and set $p=0.5$. The analytical value is a good approximation only for small values of $\tau$. The mean-field approximation becomes less accurate as more connections from past time-steps are kept in memory. Thus, the analytical estimation provides only a lower bound, which together with the solution for $\tau=1$ (dashed lines) $-$that constitutes an upper bound$-$, marks the region where spreading is possible (red regions). In other words, for a general value of $\tau$, the threshold will be lower than the analytical value computed for $\tau=1$, and larger than the corresponding value computed at $\tau$.

Overall our results highlight how the spreading of computer viruses based on social engineering is critically affected by the temporal nature of our interactions and different susceptibilities to cyber threats. Our findings show that networks’ dynamics and their interplay with the characteristics of users have to be considered in order to avoid misrepresentation of the spreading power of computer viruses in social networks. We have also quantified the extent to which the previous mismatch is important for three plausible scenarios. We, however, note that we have studied a simple network model that neglects a range of properties of real social networks such as the presence of weak and strong ties, high order correlations, and community structures. The study of the impact of these features on the unfolding of computer viruses calls for additional research.

This material is based upon work supported by, or in part by, the U. S. Army Research Laboratory and the U. S. Army Research Office under contract/grant number W911NF-18-1-0376. Y. M. acknowledges support from the Government of Aragón, Spain through grant E36-17R (FENOL) and by MINECO and FEDER funds (grant FIS2017-87519-P). The authors thanks Andrea Baronchelli and Michele Starnini for useful discussions.

I Supplementary Material

Here, we provide supplemental information about the mathematical derivation and present the sensitivity analysis of the results to the variation of the main parameters.

II Analytical derivations

We consider two types of viruses. The first does not have access to the list of contacts of each victim, thus spread only via the connections activated at each time step $t$. The second instead, is able to access the list of contacts in the last $\tau$ time-steps. As we will see in details later, at the mean-field level, the structure of the equations regulating the variation of the number of infected individuals at early times is the same. However, for simplicity, let us consider first the first type of virus.

Each node is assigned to one of $Q$ categories $x\in[1,\ldots,Q]$ that distinguish nodes according the their gullibility $\lambda_{x}$, which describes the probability that they will not recognize the threat for example clicking on the added piece of malicious content, and the time they need to recover for a successful attack, $\mu_{x}^{-1}$ ($\mu_{x}$ is the recovery rate). Nodes are also characterized by their activity $a$ which describes their propensity to engage in social interactions in the unit time. To account for observations in real systems, we extract activity from a heterogenous distribution in particular a power-law $F(a)=Ba^{-\alpha}$ with $a\in[\epsilon,1]$ to avoid divergences. Each active node creates $m$ random connections. The target of each communication act is selected with probability $p$ within nodes in the same category and with probability $1-p$ with nodes in other categories. In both cases the actual target is selected at random. The virus will be able to spread only via the connections created by active and infected users. In particular, suppose that node $i$ has been compromised. At time $t$ the node activates and sends $m$ legitimate messages to $m$ users. During the same time-step, the virus, unbeknownst to $i$, will send a message with malicious content to all $m$ users. In these settings, the variation of number of infected nodes in each activity class $a$ and category $x$ can be written as:

[TABLE]

The first term on the right hand side, describes the recovery process. The second term instead describes susceptible nodes that are connected by active nodes in the same category that are infected. These nodes get infected with probability $\lambda_{x}$ and selected with probability $p\frac{m}{N^{x}}$ (where $N^{x}$ is the total number of nodes in the category $x$). The third term, accounts for the same process but in which the susceptible node in activity class $a$ receives a message from active and infected nodes in other categories. Each node is selected with probability $(1-p)\frac{m}{N-N^{y}}$ by active vertices in class $y$. At early stages of the spreading we can assume that the number of infected to be very small respect to the susceptible thus we can approximate $S_{a}^{x}\sim N_{a}^{x}$. This is equivalent to neglect terms of the order of $(I_{a}^{x})^{2}$. We can also define $\int da^{\prime}a^{\prime}I_{a^{\prime}}^{x}=\Theta^{x}$, thus summing over all activity classes we get:

[TABLE]

In order to characterize the behavior of the number of infected at such early times, we can write, starting from Eq. 3 the equation for each auxiliary function $\Theta^{x}$. In particular, we can multiply both sides of Eq. 3 for $a$ and integrate over all classes of activity. Doing so, we obtain:

[TABLE]

where we have multiply and divided the third term for $N^{x}$. We can now define $F_{x}(a)=\frac{N_{a}^{x}}{N^{x}}$ as the distribution of activities in the category $x$, and thus $\int da\frac{aN_{a}^{x}}{N^{x}}=\int daaF_{x}(a)=\langle a\rangle_{x}$ is the average activity in the category. Finally, we let’s define $c_{x,y}=\frac{N^{x}}{N-N^{y}}$ which is acts as the mixing probability between categorie. In these settings we get:

[TABLE]

Thus we have a system of differential equations made of $2Q$ equations. In particular, we have two equations for each $x$ in the form:

[TABLE]

The conditions for the spreading can be identified by studying the eigenvalues of the Jacobian matrix of such system. The Jacobian can be written as follows:

[TABLE]

Substituting the general terms with the actual partial derivatives we get:

[TABLE]

where we defined $\beta_{x}=m\langle a\rangle_{x}\lambda_{x}$. It is important to notice the peculiarities of the Jacobian. The first $Q\times Q$ block made of the partial derivatives of the $g^{x}$ functions in the various $I^{x}$ is a diagonal block that features the recovery rates of each category. The second block on the bottom left side is a $Q\times 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 the diagonal the same function $-\mu_{x}+p\beta_{x}$. Due these properties, $Q$ eigenvalues are negative and equal to the negative of each recovery rate. The largest eigenvalue instead can be written as

[TABLE]

where $\Xi$ is an algebraic term function of all the $\beta_{x}$, $\mu_{x}$ and $c_{x,y}$. We focus on $\Lambda_{max}$ because the virus will be able to spread if and only if the largest eigenvalue is larger than zero. From this observation we obtain the conditions spreading:

[TABLE]

where $R_{0}$ is the reproductive number defined as the number of infected nodes generated by an initial seed in a fully susceptible population. It is important to mention that for any number of categories $\Xi$ has an analytical expression. However, since it derives from the characteristic equation of the Jacobian matrix, $\Xi$ gets more and more complicated as the dimensionality of the matrix increases. Generally speaking for $Q$ categories $\Xi$ is a polynomial of order $Q$ in all variables.

II.1 Q=1

In case of single category the expression of $R_{0}$ becomes:

[TABLE]

In fact, in this limit $p=1$ and the Jacobian matrix reduces to

[TABLE]

The two eigenvalues are $-\mu$ and $-\mu+\beta$. Thus the disease will be able to spread only if $\beta>\mu$.

II.2 Q=2

In the case of two categories, $Q=2$, the Jacobian becomes:

[TABLE]

In these settings we have:

[TABLE]

It is important to notice how with two categories, independently of their sizes $c_{1,2}=c_{2,1}=1$. In fact, the two sizes are constrained by $N=N^{1}+N^{2}$. Thus we have:

[TABLE]

The expression of $\Xi$ reduces to:

[TABLE]

II.3 $Q>2$

As mentioned above, in the most general case of $Q$ categories, the expression of $\Xi$, becomes quite complex. However, its expression is set unequivocally by the characteristic equation of the Jacobian matrix and can be easily obtained with any programming language that allows symbolic computations such as Mathematica. The problem can be significantly simplified in case some of the variables describing the system are set. For example in the case of $Q=3$ one might wonder what is the critical value of $\lambda_{1}$ in a system in which $\beta_{y}$ (with $y=[2,3]$) and $\mu_{y}$ with ($y=[1,2,3]$) are set. In these settings, as shown later on, it is extremely easy to compute the largest eigenvalue of the Jacobian for the particular system under consideration as function of $\lambda_{1}$.

II.4 $\tau>1$

We now turn the attention to the second type of virus that is able to access not only the connections establish at time $t$ but also those in previous $\tau$ time steps. In order to characterize the conditions for the spreading in this case, let us first understand how many people del virus will be able to reach from each node of activity $a$. This number is equal to the out-degree of those nodes. In the case considered in the previous sections $\tau=1$, thus the virus was able to reach only the nodes contacted by each active and infected node within the time-step $t$. By construction, the out-degree of such nodes is $k^{out}(a)=ma$, since their are active with probability $a$ and when active they create $m$ random connections. What about for $\tau=2$?Active nodes at time $t$ might either have $m$ connections or $2m$. The first group describes nodes that were not active at time $t-1$ but they were active at time $t$. The second group instead describe nodes that were active in both time steps. Thus:

[TABLE]

In fact, nodes of activity $a$ are not active with probability $1-a$ and are active two times in a row with probability $a^{2}$ (since the events are independent). The same reasoning applies for $\tau=3$. Here we could have three groups having either degree $m$, $2m$, and $3m$. As before, the first group describes nodes that were not active at time $t-2$ and $t-1$ but they were active at time $t$. The second group instead accounts for all the nodes that were active two times. Finally the third those that were active three times. Thus we get:

[TABLE]

In the case $\tau=4$ instead we have:

[TABLE]

It is clear that the structure of the out-degree for a general $\tau$ can be written as:

[TABLE]

Within a mean-field approximation, we can approximate the process assuming that the virus will try to infected $k^{out}(a)$ other nodes as for the case $\tau=1$. This is an approximation because each active node, at time $t$, as a quenched list of contacts, those established in the time-steps before. The node will not re-draw them ex novo as in the case $\tau=1$. Thus, we can expect the approximation to be closer to the actual process for small values of $\tau$. Within such approach, the structure of the equation is the same as those above, the only different is in the $\beta$s since we will have $m\langle a\rangle_{x}\rightarrow m\left[\langle a\rangle_{x}+(\tau-1)\langle a\rangle_{x}^{2}\right]$.

III Features of the phase space for $Q=2$

Let’s first consider the case in which $\mu_{1}=\mu_{2}=\mu$. The expression of $R_{0}$ reduces to:

[TABLE]

In the limit $p=0$, nodes in each category will connect just with nodes in the other. The expression of $R_{0}$ becomes: $R_{0}=\sqrt{\beta_{1}\beta_{2}}/\mu$. In the opposite limit, $p=1$, nodes in the two categories are separated. Thus we have two independent conditions that have the same mathematical form we encountered for $Q=1$. In fact, we have $R_{0}^{1}=\beta_{1}/\mu$ and $R_{0}^{2}=\beta_{2}/\mu$. The virus will be able to spread in the system in case either of the $R_{0}^{x}$ are larger than one. Of course, in case both are larger than one each group will experience the virus. What happens in case $0<p<1$? It is interesting to notice how the value of $R_{0}$ for a general $p$ is bounded by the $R_{0}^{x}$ of the two categories taken in isolation: $\min_{x}R_{0}^{x}\leq R_{0}(p)\leq\max_{x}R_{0}^{x}$. Before the mathematical proof, let us try to develop the intuition behind. Suppose that $\beta_{1}>\beta_{2}$. Any value of $p<1$, will reduce the spreading power of nodes in the first category. In fact, nodes in category one will be connected to some nodes in category two that are less gullible, or less active, or create a smaller number of connection (remember that $\beta_{x}=m\langle a\rangle_{x}\lambda_{x}$). Conversely, nodes in category two, will get in contact with nodes that increase the spreading potential of the virus. In order to prove this, let us consider the case $\beta_{1}>\beta_{2}$. We have to show how $R_{0}^{1}>R_{0}(p)$ and $R_{0}^{2}<R_{0}(p)$. Let us consider the first condition:

[TABLE]

which is equivalent to:

[TABLE]

This condition is respected in case $\beta_{1}(2-p)-p\beta_{2}>0$, $p^{2}(\beta_{1}+\beta_{2})^{2}+4\beta_{1}\beta_{2}(1-2p)>0$ and $(\beta_{1}(2-p)-p\beta_{2})^{2}>p^{2}(\beta_{1}+\beta_{2})^{2}+4\beta_{1}\beta_{2}(1-2p)$. The first condition implies $\beta_{1}>\frac{p\beta_{2}}{2-p}$, which is always true since $\beta_{1}>\beta_{2}$ was the initial assumption. Furthermore, it is easy to show that equation $p^{2}(\beta_{1}+\beta_{2})^{2}+4\beta_{1}\beta_{2}(1-2p)=0$ as no solution in $p$, thus the condition is always respected. Finally, the third condition implies

[TABLE]

that reduces to $\beta_{1}>\beta_{2}$. The three conditions prove Eq. 23 for all $p$. We have now to prove

[TABLE]

which is equivalent to:

[TABLE]

This condition is respected in region in which $\beta_{2}(2-p)-p\beta_{1}\geq 0$, $(\beta_{2}(2-p)-p\beta_{1})^{2}<p^{2}(\beta_{1}+\beta_{2})^{2}+4\beta_{1}\beta_{2}(1-2p)$ and $p^{2}(\beta_{1}+\beta_{2})^{2}+4\beta_{1}\beta_{2}(1-2p)\geq 0$, $\beta_{2}(2-p)-p\beta_{1}<0$. The first two conditions are respected when in the region $\frac{p\beta_{1}}{2-p}\leq\beta_{2}<\beta_{1}$. The other two instead in the region $\beta_{2}<\frac{p\beta_{1}}{2-p}$. Overall, Eq. 26 is valid in the union of these two that implies $\beta_{2}<\beta_{1}$ which is exactly the initial assumption.

Let’s consider now the general case in which also the two recovery rates are different. In the limit $p=0$, we have $R_{0}=\frac{\sqrt{(\mu_{1}-\mu_{2})^{2}-4\beta_{1}\beta_{2}}}{\mu_{1}+\mu_{2}}$. In the opposite limit instead, $p=1$, the two categories are independent thus we have two conditions as before: $R_{0}^{1}=\beta_{1}/\mu_{1}$ and $R_{0}^{2}=\beta_{2}/\mu_{2}$. It is interesting to notice how in case the two recovery rates are not the same, the phase space of the process becomes significantly more complex. In fact, differences in the rate at which nodes recovers might create interesting non-linear behaviors. In particular, consider a scenario in which the first category features a larger $\beta_{1}$ and $\mu_{1}$ respect to the second. Thus, such nodes are more prone to infection but recover faster. In case $p<1$, the coupling between the two categories might boost the spreading of the virus, since the node in category one are able to infect those in two which, although less prone to the disease stay infected for longer. For a given configuration of parameters (i.e. setting $\beta$s and $\mu$s) we can analytically determine the value of $p$ above which this phenomenon is observed. In particular, let’s assume that $\beta_{1}/\mu_{1}<\beta_{2}/\mu_{2}$. Next, we need to compute the value of $p$ (if any), for which $\beta_{2}/\mu_{2}<R_{0}(p)$. This implies:

[TABLE]

that can be written as:

[TABLE]

It is important to notice how this inequality is at the first order in $p$. Indeed, all second order terms cancel out. The value of $p$ that verifies the above inequality lays in the union of two systems of inequalities: i)$\beta_{2}(\mu_{1}+\mu_{2})-\mu_{2}p(\beta_{1}+\beta_{2})<0$ and the quantity inside the square root is larger equal than zero, ii) $\beta_{2}(\mu_{1}+\mu_{2})-\mu_{2}p(\beta_{1}+\beta_{2})>0$, and $(\beta_{2}(\mu_{1}+\mu_{2})-\mu_{2}p(\beta_{1}+\beta_{2}))^{2}<\mu_{2}^{2}\Xi^{2}$. Extensive numerical computations show that the values inside the square roots are always positive. Furthermore, the first condition in the first system result in values of $p$ always larger than one. Thus, the first system does not provide any physical ($p<1$) condition. Conversely, the first condition in the second system implies $p<1$ while the second:

[TABLE]

Thus, this is the only physical condition necessary to observe a reproductive number larger than in each category in isolation. Clearly, in the case $\beta_{2}/\mu_{2}<\beta_{1}/\mu_{1}$ the condition above becomes:

[TABLE]

In Figure 4 we verify the above condition. In particular, we set the values of $\beta_{x}$ and $\mu_{x}$ and plot $R_{0}$ from Eq. 11 as function of $p$. In particular, we consider that nodes are assigned to the categories randomly. The shaded area is the region where $\min_{x}\beta_{x}/\mu_{x}\leq R_{0}\leq\max_{x}\beta_{x}/\mu_{x}$. The vertical line show the 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 than the the $R_{0}^{x}$ of each category in isolation.

It is important to stress how the region of the phase space in which we observe this phenomenon is generally speaking quite limited. In fact, it might happen only in case the category with the larger recovery rates has also the larger gullibility. In Figure 5 we show as contour plots the region of parameters where the reproductive number of the system is larger than that correspondent value in the two categories in isolation. In particular, we set $\lambda_{1}=0.9$, $\lambda_{2}=0.8$ (A), $\lambda_{2}=0.6$ (B), $\lambda_{2}=0.4$ (C), $\lambda_{2}=0.2$ (D) and show as function of $\mu_{1}$ and $\mu_{2}$ the value of $p^{*}$. It is clear this region increases as the difference between the two gullibilities increases.

It is important to notice how the expression for $p^{*}$ is perfectly in line with the case in which $\mu_{1}=\mu_{2}$. Indeed, in this limit we get $p^{*}>1$ which implies, as expected, that the necessary condition to have a reproductive number larger than in each category in isolation is to have different recovery rates.

IV Numerical simulations

In this section we will present the sensitivity analysis to the model’s parameters. We will first consider two categories ($Q=2$). As shown in the main text, we adopted two main approaches to assign node to categories. The first is at random, the second is instead in decreasing order to activity. In particular, we order activity in decreasing order and then assign the first $gN$ nodes to the first category and the remaining to the second. Thus $\langle a\rangle_{1}=\int_{a_{c}}^{1}daaF(a)$ and $\langle a\rangle_{2}=\int_{\epsilon}^{a_{c}}daaF(a)$ and $a_{c}$ is determined in such a way that the fraction of nodes in the first class is $g$. This can be easily done imposing:

[TABLE]

Since $F(a)=\frac{1-\alpha}{1-\epsilon^{1-\alpha}}a^{-\alpha}$ we get:

[TABLE]

It is important to notice that in Eq. 11 the expression of $\langle a\rangle_{x}$ in the two assignment scenarios is slightly different. In particular we defined $\langle a\rangle_{x}=\int daF_{x}(a)a=\int da\frac{N_{a}^{x}}{N^{x}}a$. In case nodes are assigned randomly to the two categories we have that $F_{x}(a)\sim F(a)$ since $N_{a}^{x}=N_{a}/g$ and $N^{x}=N/g$ (where $g$ is the fraction of node in the general category $x$ in this case). Thus, $\langle a\rangle_{x}=\langle a\rangle$ for the two categories. In case instead nodes are assigned in decreasing order of activity $\langle a\rangle_{x}=\langle a\rangle/g$. In fact, in this limit $N_{a}^{x}=N_{a}$ (since nodes are assigned to categories as function of their activity) but $N^{x}=gN$.

In Figure 6, we consider the case of $\mu_{1}=\mu_{2}$ in case of randomly selected nodes for $p=0.9$ (A-D) and $p=0.4$ (B-E). We also consider the case of correlation between category and activity for $p=0.4$ (C-F). Respect to Figure 1 of the main text, we used a different value of $\lambda_{2}=0.3$. Across the board the analytical solution match perfectly the numerical estimation of the threshold.

In Figure 7 we consider the general case of different recovery rates. In particular we consider a different set of values respect to those in the main text. In particular, we set $\mu_{1}=10^{-2}$, $\mu_{2}=10^{-1}$, $\lambda_{2}=0.3$ and $m=4$. In panels A-B-D-E we consider random assignment of nodes to categories. In C-F we consider the correlation between activity and category. We assign to category one to most active nodes. Also, in panels B-C-E-F we considered $p=0.6$ while in panel A-D we set $p=0.9$. Overall, the figure confirms the validity of the theoretical approach and highlights one more time the effects of correlations between category assignment and activity that reduce the non-active phase space (see panel F). Furthermore, it is important to notice how the critical value in case of a single category with a recovery rate average of the two here would be $\lambda_{1}^{c}=2.5$ (not shown in the figure) which implies that the virus would not be able to spread since all the gullibilities should be smaller or equal to $1$. This confirms the importance of accounting for the presence of different categories of users in order to correctly capture the spreading power of the virus.

In Figure 8 we test the sensitivity to the parameter $m$. In the main text as well as in many of the other plots we set $m=4$. Here, we fix instead $m=6$ keeping all the other parameters the same as in the Figure 7. The analytical solutions one more time match the numerical simulations and the contour plots confirm the picture discussed in the main text and all the other similar plots.

In Figure 9 we test the sensitivity to the exponent of the activity distribution. In all the other plots we set $\alpha=-2.1$, here instead we consider $\alpha=-2.5$. We considered a scenario in which the recovery rates of the two categories is the same, set $\lambda_{2}=0.4$, $m=6$ and consider two different values of $p$. As clear from the figure, also in this case the analytical estimation matches the numerical simulations. Furthermore, it is interesting to notice how, in case of faster decay of the activity distribution (i.e. smaller value of the exponent $\alpha$), 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 gets closer to the average activity of the whole network.

IV.1 $Q=3$

Here we consider the case of three categories. For simplicity let’s consider nodes are assigned to the categories at random and that categories have the same size $N^{x}=N/3$. Also, let’s us set the values of $\beta_{x}$ with $x=[2,3]$, $\mu_{x}$ with $x=[1,2,3]$, $m=4$, and assume that links are created randomly between categories thus $p=1/3$. In particular, if we set $\beta_{2}=\beta_{3}=0.3$, $\mu_{1}=\mu_{2}=\mu_{3}=0.01$, we can use Eq. 9 to obtain the critical value of $\lambda_{1}$. In particular, the general expression of the Jacobian is:

[TABLE]

Since the categories have the same size:

[TABLE]

Plugging all the values and solving for $\lambda_{1}$ we obtain:

[TABLE]

In Figure 10 we show the comparison between the analytical prediction and the numerical simulations which perfectly matches.