Directionality reduces the impact of epidemics in multilayer networks

The study of how diseases spread has greatly benefited from advances in network modeling. Recently, a class of networks known as multilayer graphs has been shown to describe more accurately many real systems, making it possible to address more complex scenarios in epidemiology such as the interaction between different pathogens or multiple strains of the same disease. In this work, we study in depth a class of networks that have gone unnoticed up to now, despite of its relevance for spreading dynamics. Specifically, we focus on directed multilayer networks, characterized by the existence of directed links, either within the layers or across layers. Using the generating function approach and numerical simulations of a stochastic susceptible-infected-susceptible (SIS) model, we calculate the epidemic threshold for these networks for different degree distributions of the networks. Our results show that the main feature that determines the value of the epidemic threshold is the directionality of the links connecting different layers, regardless of the degree distribution chosen. Our findings are of utmost interest given the ubiquitous presence of directed multilayer networks and the widespread use of disease-like spreading processes in a broad range of phenomena such as diffusion processes in social and transportation systems.

epidemic threshold as a function of the directionality and the coupling strength between layers. In addition, we analytically derive the epidemic thresholds using generating functions, which allows to provide theoretical insights on the underlying mechanisms driving the dynamics of these systems. Our results show that the presence of directed links results in larger epidemic thresholds with respect to the case of undirected networks, and that the system is more resilient when the interlayer links are directed.

I. RESULTS
We first present results of numerical simulations of a stochastic susceptible-infected-susceptible model. In this model, nodes can be either susceptible or infected. The latters spread the disease to the formers if they are in contact with a given probability. One of the main characteristics of multiplex networks is the existence of several types of links. Thus, it is possible to associate different spreading probabilities to each of these links [14]. In our model, we assume two spreading probabilities: the interlayer spreading probability, γ, and the intralayer spreading probability, β. Hence, an infected node transmits the disease with probability β to those susceptible neighbors of the same layer and with probability γ to those located in other layers (for further details see Materials and Methods, section III B). This distinction implies that it is possible to find a critical value of β for each value of γ and vice versa. Thus, henceforth we will define the epidemic threshold as β c and explore its value as a function of γ.
The SIS dynamics is implemented on directed multiplex networks composed by two layers. As previously stated, we explore four different configurations of directionality denoted as DUD, DDD, UDU and UUU. Furthermore, to define the degree distribution in the layers we use power-law and Poisson distributions, which correspond respectively to Scale-Free (SF) and Erdős-RÃ c nyi (ER) network models. In figure 1 we show the evolution of the epidemic threshold, β c , as a function of γ for the configurations with undirected interlinks, UUU and DUD, both for ER (1A) and SF (1B) networks, for two different average degrees k .
For the cases in which the interlinks are directed, we need to define how many links point from one layer u to another layer v, either in the u → v direction or in the opposite one, u ← v. Indeed, if we set all interlinks to have the same direction, the epidemic threshold would be trivially the one of the source layer and thus the multiplex structure would play no role. For this reason, for each directed link connecting layers u and v we set the directionality to be u → v with probability p and u ← v with probability (1 − p). Consequently, in networks with directed interlinks the epidemic threshold will be given as a function of this probability p. We refer to this procedure of generating interlinks as the p-model. The same dependence of the critical threshold depicted in Figure 1 is shown in Figure 2 for DDD and UDU configurations built using the p-model.
It is also possible to study scenarios in which each interlink does not only have one possible directionality, either u → v or u ← v, but instead are bi-directional. This is achieved by setting two independent probabilities −one for each direction −, thus allowing for the coexistence of single directionality and bi-directionality in the interlinks. This situation, which we denote as the pq-model, is further analyzed in the Supplementary Information, section 1.
Lastly, in order to obtain insights into the mechanisms behind driving the spreading process on the directed multiplex networks, we analytically derive the epidemic threshold for all the configurations considered in this work, both for ER and SF networks. To this end, we extend the generating function formalism, which has been used previously in the context of directed monolayer networks [15] and interdependent directed networks [16], to multiplex networks. This formalism is outlined in the next subsection I A. The complete derivation of the epidemic threshold is presented in the Supplementary Information, section 1.

A. Generating function
Within the generating function formalism, a node has an in-degree j, out-degree k and inter-degree m with probability p jlm , being the first two related to the links contained in each layer and the latter to links connecting nodes in different layers. The generating function for the degree distribution of a node is then defined as so that in order to describe a particular network it is sufficient to set p jlm to the degree distribution of the network. Indeed, with this function it is possible to characterize several properties of the network such as the excess degree which is the main quantity that is needed for the derivation of the epidemic threshold. The excess degree of a node is defined as the number of links of a node reached by following a randomly chosen link, without including the incoming link. Hence, the distribution of excess degree of a node that is reached by following a directed link in its direction is generated by where k d is the average directed degree and the superscript (1, 0, 0) refers to partial derivation with respect to x. Similar expressions can be obtained for the excess degree of a node reached via the reverse direction of the same directed link and via an undirected link (see Supplementary Information, section 1). The size of an outbreak, as well as the epidemic threshold, can be obtained by computing the fraction of occupied links in the network. In this context, occupied link refers to a link through which the disease was transmitted. This can be accounted for by incorporating the transmissibility, i.e., the mean probability of transmission between individuals [17], to the previous equations so that where T and T uv denote the transmissibility within a layer and across layers, respectively. Recalling that β is the within-layer transmission rate, that γ is the transmission rate through links that connect different layers and that µ is the recovery rate, these quantities can be expressed as (see Supplementary Information, section 1F):  The generating function for the distribution of the size of an outbreak can be expressed as where h 1 and h 12 are recursive functions that generate the distribution of the size of an outbreak starting at a link connecting nodes in layer 1 and at a link connecting nodes in layer 1 and 2 respectively. The average size of an outbreak will be then given by the derivative with respect to w of g(w; T, T uv ) evaluated at w = 1. The said derivative goes to infinity when its denominator equals 0, which characterizes a phase transition from a phase in with only small size outbreaks to one characterized by the occurrence of macroscopic outbreaks. Thus, the epidemic threshold can be obtained from the equality and H i refers to transmission within layer i and H ij to transmission from layer i to layer j.
The above expression is general enough as to be used in the calculation of the epidemic threshold for each of the cases considered in this work. To this end, the only step that is left is to substitute H i by H d if the links in layer i are directed or, conversely, by H u if they are undirected (see Supplementary Information, section 1 and figures S1-S2 ).
In what follows, we present the results obtained for the thresholds after considering directionality (or lack thereof) and different network topologies.

B. ER networks
In ER networks the degree distribution follows a Poisson distribution. If we consider an UUU network with nodes in both layers following said degree distribution, the generating function, (1), is Inserting this expression in (7), the epidemic threshold can be expressed as (the full derivation is presented in Supplementary Information, section 1) Henceforth, to facilitate readability and unless otherwise stated, we provide expressions for the epidemic threshold in terms of the average transmission probability through intralinks, T , and the average transmission probability through interlinks, T uv . Nevertheless, the thresholds can be easily rewritten in terms of β c in a straightforward way using (4) and (5). For this case, we can rewrite the above equation and explicitly express the value of β c as, Note that if we set γ = 0 in (5) so that the spreading from one layer to the other is completely removed, T uv = 0 and (9) is simplified to βc µ = k −1 , which is the classical value of the epidemic threshold in single layer ER networks [18].
In a DUD network with nodes in both layers following a Poisson degree distribution, with the same average degree for both incoming and outgoing links, the generating function (1) is Again, by inserting this expression in (7) we obtain On the other hand, using the p-model previously described, the epidemic threshold in DDD configurations as a function of p is with m = p(1 − p)T 2 uv and in the UDU configuration is figure 3A we compare the behavior of these four configurations plotting β c as a function of γ.

C. SF networks
In SF networks the degree distribution follows a power-law of the form P (k) ∼ k −α . Thus, the epidemic thresholds are for the UUU configuration, for the DUD configuration, with m = p(1 − p)T 2 uv for the DDD configuration and with m = p(1 − p)T 2 uv for the UDU configuration. The full derivation can be found in the Supplementary Information, section 2. As in the ER case, the explicit dependence of β c with γ is shown in figure 3B.

II. DISCUSSION
Our results show that directionality is a key factor in the spreading of epidemics in multiplex networks. Even more, these findings suggest that its effects cannot be trivially generalized as the consequences of changing the directionality of some links are completely different for Scale-Free and Erdős-RÃ c nyi networks. In particular, in figure 1A, we can see that for networks with k = 6 the epidemic threshold is very similar in both UUU and DUD configurations. This effect is again seen for denser networks, k = 12, implying that it is the directionality of the interlinks, and not the one of the links contained within layers, the main driver of the epidemic in these networks. On the other hand, in figure  1B we can see that this behavior is not replicated for SF networks. Certainly, there is a large difference between the curves of the UUU and DUD configurations, implying that the directionality of intralinks is much more important in this type of networks. In agreement with these observation, when the interlinks are those that are directed, we found the same difference between ER and SF networks. As can be observed in figure 2A, the evolution of the epidemic threshold as a function of γ is again quantitatively similar for both DDD and UDU configurations. Conversely, in figure 2B, a difference between these configurations arises again for SF networks. Besides, in all the cases considered so far, figures 1 and 2, the epidemic threshold is always lower for those configurations with undirected links within the layers, compared to those in which those links are directed, given the same interlink directionality.
To get further insights into the mechanisms driving the behavior observed previously, we rely on the analytically derived thresholds and explore the evolution of β c as a function of γ for the whole range of possible values of the latter parameter. Results are shown in figure 3. In this case, we can see that the value of the epidemic threshold of the DUD configuration in SF networks tends to the value of the UUU case for large values of the spreading probability across layers, mimicking the behavior of ER networks. Thus, when γ → 1 we reach the state in which both networks exhibit the same properties, namely: (i) the epidemic threshold in DUD and UUU configurations is the same; (ii) XDX (X=U or D) configurations are almost not affected by the value of γ, except for the weakly couple regime (i.e., small values of γ). Hence, in general, one can conclude that the directionality (or lack of) of the interlinks is the main driver of the epidemic spreading process. The exception is the limit of small spreading from layer to layer, as in this scenario, the directionality of interlinks makes SF networks much more resilient, see the dashed-dotted line in 3B. Altogether, the general conclusion is that directionality reduces the impact of disease spreading in multilayer systems.
It is important to note that these results are not only relevant for the situations described in the introduction of this paper. First, because even though a system might be commonly presented as a monolayer network, it may be possible to detect different types of links in the network that would allow for the construction of a multiplex network. If this is done, as we have shown in this paper, the definition of the directionality of the interlinks is far from trivial as it can have dramatic consequences on the dynamics. In particular, the epidemic threshold can change by up to a factor of two depending on the directionality of the interlinks. Even more, these results are not restricted only to epidemic modeling, as these kind of diffusion processes can be applied to a broad range of systems. For example, the generating function approach has been proposed as a tool to identify influential spreaders in social networks [19].
One particularly interesting and open challenge is to quantify the effects that the interplay between different social networks could have on spreading dynamics. The theoretical framework developed here is particularly suitable to study this and similar challenges related to the spreading of information in social networks. On the one hand, because social relations are, by default, directed: a user is not necessarily followed by her followings, i.e., social relations are not always reciprocal [20]. On the other hand, disease-like models have been widely used to study information dissemination, or in other words, simple social contagion [21,22]. We have analyzed the dependence of the epidemic threshold with the inter-spreading rate in a real social network composed by two layers, see figure 4A. The first layer of the multilayer systems is made up by the directed set of interactions in a subset of users of the now defunct FriendFeed platform, whereas the second layer is defined by the directed set of interactions of those same users in Twitter. Even though this multiplex network originally corresponds to a DUD configuration, we have also explored the other possible configurations for the directionality of the links. Note that in contrast with the synthetic networks studied in the previous section, in this network the layers have different average degrees. In particular, the FriendFeed layer has 4,768 nodes and 29,501 directed links, resulting in an average out-degree of 6.19, and the Twitter layer is composed by 4,768 nodes and 40,168 directed links, with an average out-degree of 8.42. Nevertheless, their degree distributions are both heavy tailed, although the maximum degree in the FriendFeed network is much larger than in the Twitter network. For details on how this network was obtained, we refer the reader to the original source of the data [23].
The results, figure 4B, confirm our findings for synthetic networks. In particular, for the range of γ under consideration, the configurations with some directionality are always more resilient against the disease. These results would imply that information travels much more easily in undirected systems than in directed systems. For instance, one could build up a directed multiplex network using Instagram and Twitter data, either in a DUD configuration if it is assumed that the likelihood of someone sharing the information from one platform to the other is independent of the source or in a DDD configuration if the likelihood of sending it from Instagram to Twitter is deemed to be different than from Twitter to Instagram. On the other hand, undirected social platforms such as Facebook and Whatsapp should be modeled using UDU or UUU configurations. According to our results, information would spread more easily through these platforms, which could be worrisome as they have recently been identified as one of the main sources of misinformation spreading [24]. Lastly, it would be possible to build similar directed multiplex networks in transportation systems [25]. In these systems, the interlinks can be modeled as undirected or directed, depending on the purpose of the study. If one is interested in taking into account the fact that, for example, a metro station can be overcrowded in the incoming direction but not in the outgoing direction, such as during the morning peak time, or the other way around, during the evening peak time, it would be necessary to use directed links. On the other hand, if congestion is not relevant for the study, those links could be regarded as undirected.
In summary, we have developed a framework that allows studying disease-like processes in multilayer networks. This represents an important step towards the characterization of diffusion and spreading processes in interdependent multilevel complex systems. Our results show that directionality has a positive impact on the system's resistance to disease propagation and that the way in which interdependent (social) networks are coupled could determine their ability to spread information. Our results could be applied to a plethora of systems and show that more emphasis should be put in studying the role of interlinks in diffusion processes that take place on top of them.

A. Multilayer networks
Multilayer networks are an extension of classical contact networks in which nodes are assigned to a given layer, u, and can be connected to nodes in the same layer or in other layers. As a result, it is possible to distinguish two types of links: intralayer links, which connect pairs of nodes in the same layer, and interlayer links, which connect pairs of nodes in different layers. This formulation is used to encode features that characterize the nodes or the links that would be otherwise hidden, such as different types of interactions in protein networks or the multiple transportation modes present in mass transit systems [8]. In particular, in this work we focus on two layer directed multiplex networks. That is, networks composed by two layers in which links, either within layers or to other layers, can be directed. Furthermore, the term multiplex, in contrast to multilayer, implies that a node can only be connected to its counterpart in the other layer. In other words, it is not possible to have more than one link in each node going to the other layer [26].

B. Stochastic simulations
The SIS dynamics is implemented on two layer multiplex networks with ER and SF topologies in the layers. In the simulations, all the nodes in the system are initially susceptible. The spreading starts when one node is set to the infectious state. Then, at each time step, each infected node spreads the disease through each of its links with probability β if the link is contained in a layer and with probability γ if the link connects nodes in different layers. Besides, each infected node recovers with probability µ at each time step. The simulation runs until a stationary state for the number of infected individuals is reached.
To determine the epidemic threshold we fix the value of γ and run the simulation over multiple values of β, repeating 10 3 times the simulation for each of those values. The minimum value of β at which, on average, the number of infected individuals in the steady state is greater than one determines the value of the epidemic threshold, β c /µ. This procedure is then repeated for several values of γ to obtain the dependency of β c with the spreading across layers. Lastly, this dependency is evaluated for 10 2 realizations of each network considered in the study and their β c (γ) curves are averaged. To understand the dynamical spreading of epidemics on directed multilayer networks, we mainly investigate how the epidemic threshold is influenced by the directionality between interacting individuals. In this section, we analytically derive, following the mathematical technique of generating functions, the epidemic threshold for SIS epidemic model on directed multilayer networks.
Consider a directed multiplex network consisting of two layers interconnected by interlinks. The directed contact between an infective individual to a susceptible individual can be within the same layer, or across different layers or a mixed of both. Depending on the directionality of links within layers and the directionality of links interconnecting different layers, we analyze all possible combinations which are (i) directed layers and undirected interlinks, denoted as DUD, (ii) directed layers and directed interlinks, denoted as DDD, (iii) undirected layers and directed interlinks, denoted as UDU and (iv) undirected layers and undirected interlinks, denoted as UUU.
For a general directed multilayer network, a node has an in-degree j, out-degree k and inter-degree m with probability p jkm . The generating function for the degree distribution of a node is defined as where G(1, 1, 1) = j,k,m p jkm = 1 satisfying the probability property.
Another quantity related to the nodal degree distribution is called the excess degree distribution, which is the distribution of degrees of nodes reached by following a randomly chosen link. The probability to reach a node is biased by nodal degrees because nodes with a higher degree have a higher probability to be chosen. The probability to reach a node by following the direction of a randomly chosen link, i.e., in-link of the reached node, is jp jkm j,k,m p jkm . The corresponding generating function for the excess in-degree j − 1, out-degree k and inter-degree m reads Analogously, the generating function for a node reached by following the reverse direction of a randomly chosen directed link, i.e., out-link of the reached node, follows and similarly, the generating function for a node reached by following an undirected inter-link reads To account for the probability of a link being infected by a disease that is transmitted from an infective individual to a susceptible individual, we further modify the generating functions. Denote T i , i ∈ 1, 2 as the average probability that a susceptible individual will be infected by an infectious individual in the same layer. Denote T uv as the average probability that an infectious individual from layer u will transmit the disease to a susceptible individual in layer v. We omit the subscript of T i when there is no ambiguous. The generating function for the distribution of the number of infected links of a randomly chosen node is obtained by incorporating the probability of disease transmission in the generating function of degree distribution, which reads Analogously, we derive the generating functions, for the distribution of the number of infected links of a node reached by following a randomly chosen directed link in the designed direction, as and similarly of a node reached by following a randomly chosen undirected inter-link, as A number of nodes can be infected starting from a single infected node within the directed multilayer network. Due to the randomness of disease spreading and the variability of contacts, the size of a disease outbreak is a random variable. To eventually determine the epidemic threshold, we first investigate the distribution of the size of an outbreak starting from a single infected node and its corresponding generating function.
Denote Pr[S = s] as the probability of the size s of an outbreak starting from a single infected node. The generating function for the size distribution is defined as g(w, T, T uv ) = s Pr[S = s]w s . To solve the average size of an outbreak, we further define the generating function for the size of an outbreak starting from a node reached by a randomly chosen directed link in the designed direction, which denotes as h(w, T, T uv ) = t Pr[S = t]w t . By adding subscript u or uv to the generating function h(w, T, T uv ), we distinguish a randomly chosen link within a layer u, u = 1, 2, and a randomly chosen interlink uv connecting layers u and v.
Starting from a single infected node reached by following a randomly chosen intra-link (links within layers), the possible ways of future transmission are: the disease spreads along an intra-link in the same layer, it spreads along an inter-link to the opposite layer, it spreads along two intra-links, it spreads along one intra-link and one inter-link, etc. The transmission diagram is shown in Fig. S1(a). To account for all the transmission possibilities, we construct a recursive relation in the generating functions. Without loss of generality, we assume the disease spreading starting from an infected node in layer 1, the generating function satisfies a recursive relation Generating function for the distribution of the size of an outbreak w along a randomly chosen interlink satisfies a recursive relation Analogously, the spreading in layer 2 itself satisfies a recursive relation h 2 (w, T, T uv ) = wH 2 (1, h 2 (w, T, T uv ), h 21 (w, T, T uv ), T, T uv ) (S10) and h 21 (w, T, T uv ) = wH 21 (1, h 1 (w, T, T uv ), h 12 (w, T, T uv ), T, T uv ) (S11) The recursive relation of the generating functions is shown in Fig. S1(b). Similarly, generating function for the distribution of the size of an outbreak along a randomly chosen node in layer 1 follows g(w, T, T uv ) = wG(1, h 1 (w, T, T uv ), h 12 (w, T, T uv ), T, T uv ) (S12) The average size E[S] of an outbreak starting from a randomly chosen node thus can be calculated by Performing the derivative with respect to w on both sides of Eq. (S8)-(S12), the derivatives for generating functions h u , h uv and g read g (w, T, T uv ) = G( For w = 1, the derivatives of generating functions are simplified as from which we obtain Analogously, the average transmission probability of individuals between different layers reads, given that the spreading rate between layers is γ, Although the derivation of the epidemic threshold has been so far done for Poisson degree distributions within each layer, this analytical framework can be generalized to different degree distributions. Specifically, we further present the generalization to directed multilayer networks consisting of scale-free networks with power-law degree distributions. For a power-law degree distribution with an exponential cutoff, the degree distribution is written as  (c) undirected layers with ER distribution and directed interlinks; (d) undirected layers with SF distribution and directed interlinks. In both ER cases the average degree is 12 and in the SF cases the minimum degree is 10 and the exponent is 2.8, resulting in an average degree of 18.50.