Simulating SIR processes on networks using weighted shortest paths

Abstract

We present a framework to simulate SIR processes on networks using weighted shortest paths. Our framework maps the SIR dynamics to weights assigned to the edges of the network, which can be done for Markovian and non-Markovian processes alike. The weights represent the propagation time between the adjacent nodes for a particular realization. We simulate the dynamics by constructing an ensemble of such realizations, which can be done by using a Markov Chain Monte Carlo method or by direct sampling. The former provides a runtime advantage when realizations from all possible sources are computed as the weighted shortest paths can be re-calculated more efficiently. We apply our framework to three empirical networks and analyze the expected propagation time between all pairs of nodes. Furthermore, we have employed our framework to perform efficient source detection and to improve strategies for time-critical vaccination.

Figure 1

The mapping of spreading dynamics to an ensemble of weighted networks. The weights represent the propagation time delays. Each weighted network represents one realization of the stochastic dynamics from an arbitrary source node. The propagation time from the source node to any other node is equal to the shortest path length in weighted network (color coded). The dotted line represents edges with infinite edge weight.

Figure 2

Markov Chain Monte Carlo sampling. The transitions W_i,j = w(G_i → G_j) between different instances of weighted networks in the ensemble. Randomization of the weight of an edge (denoted by red) corresponds to the transition between the states (realizations of the process).

Figure 3

Toy network model example. (a) Sketch of simple toy network, that consists of n_c chains with l nodes. The probability that source node s infects destination node d can be calculated analytically (see supplementary information, Sec. 5. for details) and is given by $P(s\to d)=1-{\sum }_{j=0}^{n_c}{p}_{n,j}\cdot {\left(1-p_{\mathrm{1,1}}^l\right)}^j$. The term p_n,k, denotes transmissibility for the first neighborhood Eq. 10 for Poissionan SIR process and Eq. 9 for generalized non-Markovian process. (b) Simulation results on the toy network with n = 20, l = 3 with Poisson SIR (β = 1, γ = 1). The blue curve represents the estimations with the exact mapping (conditional independence among edge weights), the red curve represents the mean-field mapping (weights are independent) and the black dotted line represents the analytical solution (see supplementary information, Sec. 5. for a million size toy network experiment).

Figure 4

Equivalence to bond percolation. Average outbreak size in time estimated with the proposed mapping and comparison with bond percolation late-time outbreak size (dotted lines), see Eq. 4. In a limit of infinite time the mapping becomes equivalent to the bond percolation. Results were obtained for the SIR dynamics (continuous time) with different transmission rates β (0.3, 0.03, 0.003, 0.0003) and recovery rate γ = 0.001 and with n = 10⁴ simulations on the 4-connected two dimensional regular lattice (11 × 11).

Figure 5

The expected propagation time D_ij on empirical networks. The average propagation time is calculated as the average shortest path distances in the ensemble of weighted networks (mapped dynamics). (a) Nodes in the email network network are placed on a circle, where the angular coordinate represents the average propagation time from the source node (black) with index 500, increasing in the clockwise direction (from blue to red). The average propagation time D_ij for: (b) the email network (β = 0.01), (c) the airport network (β estimated from flux data see Supplementary information, Sec. 7), (d) the Petster network (β = 0.01). Nodes are ordered by their characteristic spreading timescale τ_i, see Eq. 8) with condition $\overline{N}=N\mathrm{/2}$.

Figure 6

Source detection. We compare source detection performance for different methods for the SIR model: (a) β = 0.7, γ = 0.3 and (b) β = 0.7, γ = 0.7 at time T = 5 (discrete time) on the 4-connected lattice (30 × 30 nodes). For each observed realization, we rank nodes according to their estimated probability score according to the Direct Monte Carlo method, in case of a tie we prioritize the node with a higher estimator according to the temporal distance method. We compare different methods as explained in the text: Direct Monte Carlo source probability ranking, source detection ranking by topological distance, and the estimation from our mapping (temporal distance). Results are averaged over 30 observed realizations and error bars show one standard deviation, from top to bottom.

Figure 7

Time critical vaccination. (a) Shows an illustration of the time critical vaccination process. Vaccination candidates are highlighted by the diamond symbols. Red note denote infected individuals, and blue nodes are susceptible. In (b) and (c), we show the total number of infected nodes up to time t for the different vaccination strategies described in the text (m = 0.2 of the total number of nodes). The first dashed vertical line denotes T₀ = 3, and the second one T₀ + τ = 13. In (b) for the discrete time SIR β = 0.03,γ = 0.01 and in (c) for the discrete time SIR β = 0.05,γ = 0.01 with source node 10 in the Petster network.