SIR Epidemic Simulator on Multiplex Graphs

This module provides efficient, reproducible SIR (Susceptible-Infected-Recovered) epidemic simulators for multiplex networks (same node set, multiple edge layers).

Two APIs Available:

OOP-Style API (recommended): Clean object-oriented interface with SIRDynamics class
Function-Based API: Lower-level functions with sparse matrix input

Features

Discrete-time simulation: Step-based updates with configurable time step
Continuous-time simulation: Gillespie algorithm for exact event-driven dynamics
Multiplex support: Multiple interaction layers with per-layer transmission rates
Flexible parameterization: Per-layer transmission rates, per-node recovery rates, layer weights
Exogenous importations: Optional external infection sources
Rich outputs: Incidence curves, event logs, per-layer attribution
Sparse linear algebra: Efficient handling of large networks
Reproducible: Seeded random number generation

OOP-Style API (Recommended)

The recommended way to run SIR simulations is using the SIRDynamics class, which provides a clean, book-aligned API.

Basic Example

from py3plex.dynamics import SIRDynamics
import networkx as nx

# Create a network
G = nx.karate_club_graph()

# Create SIR dynamics
sir = SIRDynamics(
    G,
    beta=0.3,              # Infection probability
    gamma=0.1,             # Recovery probability
    initial_infected=0.1   # 10% initially infected
)

# Set seed for reproducibility
sir.set_seed(42)

# Run simulation
results = sir.run(steps=100)

# Extract measures
prevalence = results.get_measure("prevalence")
state_counts = results.get_measure("state_counts")

print(f"Peak prevalence: {prevalence.max():.2%}")
print(f"Final recovered: {state_counts['R'][-1]}")

Multilayer Example

from py3plex.dynamics import SIRDynamics
from py3plex.core import multinet

# Create multilayer network
network = multinet.multi_layer_network(directed=False)

# Add nodes to two layers
for i in range(20):
    network.add_nodes([
        {'source': i, 'type': 'physical'},
        {'source': i, 'type': 'digital'}
    ])

# Add edges (simplified)
for i in range(19):
    network.add_edges([{
        'source': i, 'target': i+1,
        'source_type': 'physical', 'target_type': 'physical'
    }])

# Run SIR dynamics
sir = SIRDynamics(network, beta=0.3, gamma=0.1)
sir.set_seed(42)
results = sir.run(steps=100)

Key Methods

set_seed(seed): Set random seed for reproducibility
run(steps): Run simulation for specified number of steps, returns DynamicsResult
results.get_measure(name): Extract measures from results
- "prevalence": Fraction infected over time (array)
- "state_counts": Counts of S, I, R over time (dict of arrays)
- "trajectory": Full state history (list of dicts)
results.to_pandas(): Convert to pandas DataFrame for analysis

Model Definition

States

Each node has a single health state shared across all layers:

S (0): Susceptible
I (1): Infected
R (2): Recovered

Transmission

Infections occur along layer-specific edges with rate/probability β[α]
Edge weights can represent contact multiplicity or strength
Hazards from all layers combine:
- Continuous-time (Gillespie): Additive combination
- Discrete-time: Independent-trial complement (1 - ∏ exp(-λ_α dt))

Recovery

Infected nodes recover at rate γ (scalar or per-node)
Recovery is independent of network structure

Optional Features

Layer weights: Scale transmission rates by layer importance
Time-varying parameters: β(α, t) and γ(i, t) as callables
Exogenous infections: Poisson import rate η(t)

Installation

The SIR multiplex module is part of py3plex. Ensure you have the required dependencies:

pip install numpy scipy

Or install py3plex with all dependencies:

pip install py3plex

Function-Based API

For advanced users who need direct access to sparse matrices and detailed control, the function-based API is also available.

Quickstart

Discrete-Time Simulation

import numpy as np
import scipy.sparse
from py3plex.algorithms.sir_multiplex import simulate_sir_multiplex_discrete

# Create a simple ring network
N = 20
edges = [(i, (i+1) % N) for i in range(N)]
row, col = zip(*edges)
A = scipy.sparse.csr_matrix((np.ones(len(edges)), (row, col)), shape=(N, N))

# Run simulation
result = simulate_sir_multiplex_discrete(
    A_layers=[A],
    beta=0.5,          # Transmission rate
    gamma=0.2,         # Recovery rate
    dt=0.5,            # Time step
    steps=100,         # Number of steps
    rng_seed=42        # For reproducibility
)

# Access results
print(f"Peak prevalence: {result.I.max()} at time {result.times[result.I.argmax()]}")
print(f"Final attack rate: {result.R[-1] / N:.2%}")

Continuous-Time (Gillespie) Simulation

from py3plex.algorithms.sir_multiplex import simulate_sir_multiplex_gillespie

result = simulate_sir_multiplex_gillespie(
    A_layers=[A],
    beta=0.5,
    gamma=0.2,
    t_max=100.0,       # Maximum time
    rng_seed=42,
    return_event_log=True  # Get detailed event log
)

# Examine events
for t, event_type, node_id, layer in result.events[:10]:
    print(f"t={t:.2f}: {event_type} at node {node_id}")

Multiplex Network Example

# Create two-layer network: physical contact + online interaction
N = 50

# Layer 1: Physical proximity (localized)
edges_physical = [(i, (i+1) % N) for i in range(N)]
row, col = zip(*edges_physical)
A_physical = scipy.sparse.csr_matrix((np.ones(len(edges_physical)), (row, col)), shape=(N, N))

# Layer 2: Online network (more random)
np.random.seed(42)
edges_online = [(i, j) for i in range(N) for j in range(N)
                if i != j and np.random.random() < 0.1]
row, col = zip(*edges_online) if edges_online else ([], [])
data = np.ones(len(edges_online)) if edges_online else []
A_online = scipy.sparse.csr_matrix((data, (row, col)), shape=(N, N))

# Simulate with different transmission rates per layer
result = simulate_sir_multiplex_discrete(
    A_layers=[A_physical, A_online],
    beta=np.array([0.3, 0.1]),  # Higher transmission in physical layer
    gamma=0.15,
    layer_weights=np.array([1.0, 0.5]),  # Weight online layer less
    dt=0.5,
    steps=200,
    rng_seed=123,
    return_layer_incidence=True  # Track which layer caused infections
)

# Analyze layer contributions
from py3plex.algorithms.sir_multiplex import summarize
summary = summarize(result)
print(f"Layer contributions: {summary['layer_contributions']}")

API Reference

simulate_sir_multiplex_discrete

Discrete-time SIR simulation with synchronous updates.

Parameters:

A_layers (list[csr_matrix]): List of L sparse adjacency matrices (N×N)
beta (float | ndarray): Transmission rate (scalar or length L array)
gamma (float | ndarray): Recovery rate (scalar or length N array)
layer_weights (ndarray, optional): Layer importance weights (length L)
dt (float): Time step size (default: 1.0)
steps (int): Number of simulation steps (default: 100)
initial_state (ndarray, optional): Initial state array with values {0,1,2}
initial_infected (ndarray, optional): Boolean mask for initially infected nodes
import_rate (float | callable): Exogenous infection rate (default: 0.0)
rng_seed (int): Random seed (default: 0)
return_event_log (bool): Return event log (default: False)
return_layer_incidence (bool): Return per-layer infections (default: False)

Returns: EpidemicResult object

simulate_sir_multiplex_gillespie

Continuous-time SIR simulation using Gillespie algorithm.

Parameters: Same as discrete, except:

t_max (float): Maximum simulation time instead of dt and steps
return_event_log (bool): Default is True

Returns: EpidemicResult object

EpidemicResult

Dataclass containing simulation results:

times (ndarray): Time points
S, I, R (ndarray): Counts of each state over time
states (ndarray, optional): Full state history (discrete only)
incidence_by_layer (ndarray, optional): Infections per layer over time
events (list, optional): Event log with tuples (time, type, node, layer)
meta (dict): Simulation metadata

basic_reproduction_number

Compute approximate R₀ using spectral radius.

R0 = basic_reproduction_number(A_layers, beta, gamma, layer_weights=None)

For homogeneous networks with uniform recovery rate γ:

\[R_0 \approx \rho\left(\sum_\alpha \frac{\beta_\alpha w_\alpha}{\gamma} A^{(\alpha)}\right)\]

summarize

Extract summary statistics from simulation results.

summary = summarize(result)

Returns dictionary with:

peak_prevalence: Maximum infected count
peak_time: Time of peak
attack_rate: Final proportion recovered
time_to_extinction: When infections reached zero
total_infections: Total recovered at end
layer_contributions: Proportion of infections per layer
duration: Total simulation time

Update Rules

Discrete-Time

For each time step t → t + dt:

Infection probability for susceptible node i:

\[ \begin{align}\begin{aligned}\lambda_{i,\alpha}(t) = \beta_\alpha w_\alpha \sum_j A^{(\alpha)}_{ji} \mathbb{1}\{X_j = I\}\\p_i(t) = 1 - \prod_\alpha \exp(-\lambda_{i,\alpha}(t) \, dt)\end{aligned}\end{align} \]
Recovery probability for infected node i:

\[p_{\text{rec}}(i) = 1 - \exp(-\gamma_i \, dt)\]
Synchronous update: Compute all probabilities from state at t, then apply changes

Continuous-Time (Gillespie)

Compute total event rate:

\[\Lambda = \sum_{i \in S} \sum_\alpha \lambda_{i,\alpha} + \sum_{i \in I} \gamma_i + \eta|S|\]
Sample time to next event: \(\Delta t \sim \text{Exp}(\Lambda)\)
Select event type proportionally to rates
Update state and incrementally adjust affected rates
Record event and repeat until \(t > t_{\text{max}}\) or no infected remain

Units and Conventions

Time: Arbitrary units (days, hours, etc.). Ensure β, γ, dt/t_max are in same units.
Rates: Events per time unit. β = 0.5 per day means average 0.5 transmissions per infected-susceptible contact per day.
Probabilities: For discrete-time, computed as 1 - exp(-rate × dt)
Edge weights: Scale transmission. Weight 2.0 doubles effective contact rate.

Performance Tips

Use sparse matrices (CSR format) for all adjacency matrices
For large networks (N > 10⁴), use smaller time steps or t_max
Disable event logs for very long simulations to save memory
Vectorized operations are used internally for efficiency

Examples

See examples/advanced/example_sir_multiplex.py for additional examples including:

Comparison of discrete vs continuous time
Effect of network structure on epidemic spread
Multi-layer network analysis
Parameter sensitivity studies

References

Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry, 81(25), 2340-2361.
Keeling, M. J., & Rohani, P. (2008). Modeling infectious diseases in humans and animals. Princeton University Press.
De Domenico, M., et al. (2013). Mathematical formulation of multilayer networks. Physical Review X, 3(4), 041022.

Citation

If you use this simulator in your research, please cite py3plex:

@software{py3plex,
  author = {Škrlj, Blaž},
  title = {Py3plex: A library for analysis and visualization of heterogeneous networks},
  url = {https://github.com/SkBlaz/py3plex},
  year = {2023}
}