Threshold driven contagion on weighted networks

Weighted networks capture the structure of complex systems where interaction strength is meaningful. This information is essential to a large number of processes, such as threshold dynamics, where link weights reflect the amount of influence that neighbours have in determining a node's behaviour. Despite describing numerous cascading phenomena, such as neural firing or social contagion, the modelling of threshold dynamics on weighted networks has been largely overlooked. We fill this gap by studying a dynamical threshold model over synthetic and real weighted networks with numerical and analytical tools. We show that the time of cascade emergence depends non-monotonously on weight heterogeneities, which accelerate or decelerate the dynamics, and lead to non-trivial parameter spaces for various networks and weight distributions. Our methodology applies to arbitrary binary state processes and link properties, and may prove instrumental in understanding the role of edge heterogeneities in various natural and social phenomena.


S1 Approximate master equations on weighted networks
In this section we justify and outline the derivation of approximate master equations (AMEs) for stochastic binary-state dynamics on weighted networks. We begin with a general derivation, and later show how this framework simplifies for monotone dynamics. We identify edge types by the value of their weights, however the formalism remains unchanged by distinguishing edge types with other link properties, such as direction or color. The following derivation builds upon a formalism developed by Gleeson [1][2][3][4].

S1.1 Binary-state dynamics
The theoretical framework of the AMEs applies to stochastic binary-state dynamics on random networks, where nodes have degree k with distribution P (k). The network is assumed to be infinite and maximally random, meaning there is no correlation between k and any other graph property. We attribute to each node in the network one of two possible states, susceptible (S) or infected (I). We denote by m the number of infected neighbours of a node, with 0 ≤ m ≤ k. We refer to m interchangeably as the infected neighbour count or partial degree. A susceptible node of degree k that has m infected neighbours belongs to the set S k,m , and an infected node to the set I k,m . As such the network can be partitioned into a finite number of sets, assuming P (k) is bounded by a minimum and maximum degree k min ≤ k ≤ k max .
We may further partition S k,m and I k,m by assuming a finite number of edge types within the network, distinguished by their weight. If edge weights take one of n distinct values, it is instructive to introduce a weight vector w = (w 1 , . . . , w n ) T to store the n values w j . This could be generalized to any n-dimensional edge property vector. Further, we define the degree vector k = (k 1 , . . . , k n ) T and partial degree vector m = (m 1 , . . . , m n ) T . Here, k j and m j respectively denote the number of neighbours and the number of infected neighbours of a node that are connected by an edge of weight w j . These quantities are connected to the degree and partial degree via k = j k j and m = j m j . Moreover, S k,m and I k,m denote the set of susceptible and infected nodes, respectively, that have degree vector k and partial degree vector m.
Every node in S k,m belongs to a corresponding set S k,m , as is the case for I k,m and I k,m . Finally, the size of these sets is quantified through s k,m (t) and i k,m (t), the fraction of nodes with degree vector k that are susceptible or infected at time t, and have partial degree vector m. These quantities enumerate all possible node configurations over the course of any binary-state process. In other words, the sets S k,m and I k,m cannot be further partitioned, making s k,m (t) and i k,m (t) ideal functions of a rate equation formalism.
A dynamical process may be induced on such a network by assigning an initial state to each node, and allowing this state to evolve according to rates F k,m and R k,m per infinitesimal time step dt. The former is the rate of infection of susceptible nodes in S k,m over a time interval dt, the latter the rate of recovery from the infected state for nodes in I k,m over dt. This means all nodes sharing a degree vector k and partial degree vector m are equivalent in their rates of infection and recovery. The rate equations governing the sets S k,m and I k,m in a dynamics allowing both infection and recovery are Figure S1: Representation of the AME system in Eqs. (S1)-(S2). The system constitutes an initial value problem, solved over a set of rate equations in s k,m and i k,m , whose gain and loss terms, with their associated rates, are illustrated in the figure. The index 1 ≤ j ≤ n enumerates the gain or loss type, corresponding to the 4n + 2 ways in which a node may enter and leave class S k,m or I k,m . This is possible through the infection and recovery of a node's neighbours (vertical movements) and through infection and recovery of a node itself (horizontal movements). and where e j is the j-th basis vector of dimension n, and β s j , β i j , γ s j and γ i j the probabilities of a j-type neighbour becoming infected or recovering over an interval dt, calculated using the full system of s k,m and i k,m values.
Explicitly, the β j terms quantify the rate of infection of j-type neighbours for both susceptible and infected nodes, and while the γ j terms give the rate of recovery of j-type neighbours for both susceptible and infected nodes, and where we sum over k min ≤ k ≤ k max , all k such that j k j = k, and all m such that 0 ≤ m j ≤ k j . The values s k,m (t) and i k,m (t) combined with the degree and degree vector distributions P (k) and P (k) give us the density of infected nodes ρ(t), The initial conditions are prescribed by i k,m (0) and s k,m (0), subject to the normalisation condition As such, we have defined a closed system of deterministic equations that can be solved numerically using standard methods (Fig. S1).

S1.2 Monotone dynamics
The above derivation assumes generic infection and recovery rates F k,m and R k,m . In this section, we illustrate a solution of the AMEs particular to monotone dynamics with the example of a threshold rule for complex contagion. This may be generalized to other monotone or non-recovery dynamics, where R k,m = 0.
In dynamical processes on weighted networks, we are typically interested in the node properties relating to the edge weight. We define the strength of a node as the sum of edge weights across all neighbours, q k = k·w.
Similarly, we define the partial strength as the sum of edge weights across all infected neighbours, q m = m·w, with 0 ≤ q m ≤ q k . The infection rate for complex contagion can thus be expressed as with rate of recovery R k,m = 0. Here, p is the rate of spontaneous infection, whereby a susceptible node may become infected independently of the state of its neighbours. The threshold φ is the fraction of a node's Figure S2: Representation of the AMEs for monotone dynamics. A recovery rate R k,m = 0 implies γ i j = γ s j = 0, so the rate equations for S k,m and I k,m are characterized by only 2n + 1 gain and loss terms, in contrast to Fig. S1. total strength that must be met by the partial strength for that node to undergo induced infection. In other words, it is the fraction of a node's total received influence that must come from infected neighbours before that node itself becomes infected. Correspondingly, the master equations become (Fig. S2) The AME system (S8) is decoupled, so we may only consider the equation for s k,m when, for example, reducing the AMEs to a lower-dimensional system.

S1.3 Monotone dynamics for bimodal weight distribution
Here we analyse the simple case of a network with arbitrary degree distribution P (k) and n = 2 edge weights, w 1 , w 2 > 0, which may or may not appear with equal probability in the network. The probability distribution P (w) of a randomly chosen weight w is and 0 elsewhere, with δ ∈ (0, 1). Assuming w 1 ≥ w 2 for the sake of simplicity, δ is the fraction of strong edges in the network, and thus contributes to skewness in the weight distribution. The weight average and standard deviation are given by We may invert the linear system in Eq. (S10) to obtain the strong and weak weights, w 1 and w 2 , in terms of µ and σ, where σ ≥ 0 and µ > σ δ/(1 − δ) 1 . Thus, we can take µ and σ as parameters, and use Eq. (S11) to obtain values for the weights in the network.
As for the AME formalism in the case of a bimodal weight distribution, the weight, degree and partial degree vectors are w = (w 1 , w 2 ) T , k = (k 1 , k 2 ) T , and m = (m 1 , m 2 ) T , respectively, subject to the constraints Further, the sum over degrees, degree vectors and partial degree vectors can be written explicitly as . (S14) With Eqs. (S13)-(S14) and a given degree distribution P (k), we may write explicitly the full and reduced AME systems, solve them numerically, and explore the behaviour of the fraction of infected nodes ρ(t) as a function of all parameters.
The bimodal case is ideal as a means of illustrating how a given node may occupy a series of (k, m) classes over the course of a dynamical process. By taking the example of a node in class k = (2, 2), m = (1, 1) adhering to the infection rate F k,m (Fig. S3), we illustrate the interdependencies of various node classes and possible flows between them. This class corresponds to nodes with degree k = 4, consisting of two strong and two weak neighbours, one of each being infected. It follows that two ways in which a node may leave this class is by an additional neighbour of either type becoming infected. Similarly, two ways in which a node may enter the class is by having only one infected neighbour of either edge type, and gaining an infected neighbour of the opposite type (Fig. S4). We note that although the degree vector k of a node is fixed throughout the dynamical process, its partial degree vector m is free to change according to the number and edge-types of its infected neighbours. This allows us to attribute relative sizes to each of the (k, m) classes. A node's class is a dynamic quantity, and it is the flow of nodes through each class that we use to characterize the state of the system through s k,m (t) and i k,m (t) over time.  Figure S4: Possible (k, m) classes for k = 3, 4 with n = 2, and flows between them. Note that it is impossible for a node to move between classes of different k, since the degree vector of a node is fixed in time. In the non-recovery model of the figure, it is also impossible to make a downward transition to a class with lower m.

S1.4 Reduced AMEs
To reduce the dimension of the weighted AMEs for monotone dynamics [Eq. (S8)] in the case of a stepwise infection rate F k,m [Eq. (S7)], we need to consider system-wide quantities that are more aggregated than s k,m . We take the probability ρ(t) that a randomly chosen node is infected, i.e. the fraction of infected nodes in the network, and the probability ν j (t) that a randomly chosen neighbour (across a j-type edge) of a susceptible node is infected (see Methods). We start by proposing an exact solution for the AME system in terms of the ansatz where B kj ,mj = kj mj ρ mj (1 − ρ) kj −mj is the binomial distribution. The meaning of the ansatz in Eq. (S15) is quite intuitive and takes into account two processes. First, a susceptible node with k j edges of type j, is connected to m j infected nodes with the binomially distributed probability B kj ,mj (ν j ). Second, for q m < φq k a susceptible node does not fulfill the threshold rule and can only become infected spontaneously with probability e −pt , since the system is progressively being filled due to spontaneous infection. Considering these processes as independent leads to the product in Eq. (S15). The next step is to insert the ansatz (S15) into the AME system (S8) and derive a set of ordinary differential equations (ODEs) for the aggregated quantities ρ and ν j . Taking the time derivativeṡ k,m of Eq. (S15) (i.e. the left-hand side of the AME system) we geṫ Then, we use the infection rate of weighted contagion for q m < φq k , the ansatz (S15) and the binomial identity in the right-hand side of the AME system to obtain Equating Eqs. (S16)-(S18) as in the AME system, and separating terms for a given value of j from the rest Since the left-hand side of Eq. (S19) depends on the function ν j and its derivative only, while the right-hand side depends on the rest of the functions ν i , both sides must be equal to some constant c j . For the left-hand side, this means thaṫ For the ODE (S20) to hold regardless of the values of m j and k j , we need c j = 0. Then, the condition on ν j such that the ansatz (S15) is a solution of the AME system iṡ This ODE has the initial condition ν j (0) = ρ(0) = 0, obtained by evaluating Eq. (S15) at t = 0 and comparing with the expression B kj ,mj (0), which corresponds to an infinitesimally small initial infection randomly distributed in the network (see Methods).
The next step is to extend a general result derived by Gleeson in [2] [Eqs. (F6)-(F10) therein] to the case of weighted networks. We start by multiplying the AME system (S8) by P (k)P (k)(k j − m j ) and summing over k, k, and m, From the definition of β s j in Eq. (S3), the first term on the right hand side of Eq. (S22) may be written as As for the second term on the right hand side, when i = j the term telescopes to Eq. (S23), and for i = j it telescopes to 0. Overall, we can rearrange Eq. (S22) and obtain with d j a constant that can be determined from initial conditions. Assuming an infinitesimally small fraction of infected nodes randomly distributed in the network (see Methods), and since ν j (0) = ρ(0) = 0 and B ki,mi (0) = δ mi,0 with δ ij the Kronecker delta, we have where z j is the average number of j-type edges a node has in the network, or average j-degree. Thus, The next step is to use Eq. (S27) to find a new expression for β s j and thus write the ODE (S21) explicitly in terms of ν j . Noting that the left-hand side of Eq. (S27) is the denominator in the definition of β s j , we get where the sums qm<φq k and qm≥φq k run over all partial degree vectors m that satisfy their respective inequalities, and we have also inserted the ansatz (S15) and the binomial identity (k j − m j )B kj ,mj (ν j ) = k j (1 − ν j )B kj −1,mj (ν j ) to simplify the expression of β s j . Moreover, we may introduce the response function of the monotone, threshold-driven dynamics of our model, with f (0, 0) = 0 (a function that activates when a (k, m)-class node fulfils the threshold condition and gets infected), in order to invert the restricted sum of Eq. (S28), Overall, comparing Eqs. (S21)-(S28), we can write an explicit ODE for ν j , with ν = (ν 1 , . . . , ν n ) T , j = 1, . . . , n, and the function g j (ν, t) given by where we have defined f t = 1 − (1 − p)e −pt .
Even though Eq. (S31) is closed and in this sense equivalent to the AME system (S8), we may also derive a corresponding ODE for ρ, since we are mainly interested in the temporal evolution of the fraction of infected nodes in the network. From the definition of ρ and the AME system we havė where the second term in the right-hand side telescopes to zero. Then, we use an algebraic manipulation similar to that of Eq. (S28) to obtaiṅ Thus, the ODE for ρ is where the function h(ν, t) is given by Combining all of these results, the AME system (S8) is reduced to a closed system of n + 1 coupled, non-linear ODEs, with the quantities g j (ν, t) and h(ν, t) given explicitly by Eqs. (S32)-(S36).

S2 Combinatorial solution of parameter space boundaries
The dynamics of threshold driven contagion on weighted networks depends on the stepwise infection rate  the relative time of cascade emergence (Fig. S5). In other words, boundaries for t r in (σ, φ)-parameter space separate network configurations where the corresponding (k, m) class does and does not satisfy the threshold rule q m ≥ φq k , thus promoting or hindering spreading. In Fig. S5 we enumerate all possible boundaries for up to two infected neighbours in the case of a k-regular random network (k = 7) and a bimodal weight distribution (n = 2). Fig. S5a shows the case where one strong infected neighbour, m = (1, 0), is sufficient to cause infection for nodes with k 1 = 1, . . . , k strong neighbours. These curves are shown in red, since the corresponding node classes induce a faster cascade of spreading compared to the same process carried out on an unweighted network. Since the weight vector is w = (µ + σ, µ − σ) T for weight mean µ and skewness δ = 0.5, boundaries can be written explicitly as

S3 Comparison of numerical experiment and AME solutions
As discussed in the main text, the behaviour of threshold driven contagion over weighted networks, in Monte Carlo simulations or as predicted by numerically computing AME solutions, is remarkably consistent. We where N σ and N φ are the number of points considered in each dimension of the parameter space, and t sim r (σ, φ) and t theo r (σ, φ) are the relative times of cascade emergence for a given (σ, φ) point, measured by numerical simulations or AME solutions. This quantity is very small, M D = 2.8×10 −7 , which indicates that despite making several simplifying assumptions during the derivation of the full and reduced AME systems [Eqs. (S8)-(S37)], they provide an extremely good approximation of the spreading process with differences only due to small statistical fluctuations in finite-size numerical simulations.

S4 Other heterogeneous synthetic and real networks
In the main text we have seen that, even for heterogeneous synthetic and real world networks, threshold driven contagion strongly depends on link weights via simple mechanisms that can be understood by master equations or combinatorial arguments, and develops spreading cascades that are either faster of slower than their counterparts in unweighted contagion, depending on the values of σ and φ. Here we further support this argument by exploring two additional examples of synthetic and empirical networks. The synthetic structure is a configuration-model random network with average degree z = 7, and bimodal weight distribution with average µ = 1 and skewness δ = 0.5 (Fig. S7a). The (σ, φ)-parameter space for t r is qualitatively very