生物医学网络科学Biomedical Network Science
Example08 / 14

集合种群模型与有效距离

所属专题
网络上的传播与免疫
下载
Notebook · .py

在 M 个城市组成的随机流动网络上运行随机集合种群 SIR（城市内均匀混合，城市间按旅行概率移动感染者）。 用 Brockmann–Helbing 有效距离 d_mn = 1 - ln P_mn（P_mn 为离开 n 的旅客去往 m 的比例）求从源城市出发的最短有效路径， 比较“到达时间”与有效距离、跳数距离的相关性。 解读：流量高度不均时，跳数距离几乎相同的城市到达时间可以差很多；到达时间与有效距离的相关明显高于 与跳数距离的相关，这就是 Brockmann 与 Helbing（2013）所说的“隐藏几何”。

Output输出图

集合种群模型与有效距离：输出图 1

Output文字输出

被感染城市 149 / 149
到达时间 vs 有效距离 r = 0.88
到达时间 vs 跳数距离 r = 0.61

Code代码

# 在 M 个城市组成的随机流动网络上运行随机集合种群 SIR（城市内均匀混合，城市间按旅行概率移动感染者）。
# 用 Brockmann–Helbing 有效距离 d_mn = 1 - ln P_mn（P_mn 为离开 n 的旅客去往 m 的比例）求从源城市出发的最短有效路径，
# 比较“到达时间”与有效距离、跳数距离的相关性。
# 解读：流量高度不均时，跳数距离几乎相同的城市到达时间可以差很多；到达时间与有效距离的相关明显高于
#       与跳数距离的相关，这就是 Brockmann 与 Helbing（2013）所说的“隐藏几何”。
import numpy as np
import networkx as nx
from scipy.stats import pearsonr, spearmanr

rng = np.random.default_rng(4)
M = 150
G = nx.barabasi_albert_graph(M, 6, seed=4)
pop = (rng.pareto(1.5, M) + 1) * 5e4
W = np.zeros((M, M))
for u, v in G.edges():                       # 流量 ∝ (P_u P_v)^0.5（重力模型的简化）
    W[u, v] = W[v, u] = np.sqrt(pop[u] * pop[v]) * rng.lognormal(0, 2.0)
travel_rate = 0.002 * W.sum(1) / W.sum(1).mean()   # 每天离开城市的人口比例
Pmat = W / W.sum(1, keepdims=True)                  # P[n, m]：离开 n 去 m 的比例

H = nx.DiGraph()
for n in range(M):
    for mm in np.flatnonzero(Pmat[n]):
        H.add_edge(n, mm, d=1 - np.log(Pmat[n, mm]))
src = int(np.argmax(pop))
d_eff = nx.single_source_dijkstra_path_length(H, src, weight="d")
d_hop = nx.single_source_shortest_path_length(G, src)

beta, gamma = 0.5, 0.25
S = pop.copy(); I = np.zeros(M); I[src] = 10; S[src] -= 10
arrival = np.full(M, np.nan); arrival[src] = 0
for day in range(1, 400):
    newinf = rng.binomial(S.astype(int), 1 - np.exp(-beta * I / pop))
    rec = rng.binomial(I.astype(int), gamma)
    S -= newinf; I += newinf - rec
    travelers = rng.binomial(I.astype(int), np.minimum(travel_rate, 1))
    for n in np.flatnonzero(travelers):
        dest = rng.multinomial(travelers[n], Pmat[n])
        I[n] -= travelers[n]; I += dest
    hit = np.isnan(arrival) & (I >= 1)
    arrival[hit] = day
ok = ~np.isnan(arrival) & (np.arange(M) != src)
nodes = np.flatnonzero(ok)
print("被感染城市 %d / %d" % (ok.sum(), M - 1))
print("到达时间 vs 有效距离 r = %.2f" % pearsonr([d_eff[n] for n in nodes], arrival[nodes])[0])
print("到达时间 vs 跳数距离 r = %.2f" % pearsonr([d_hop[n] for n in nodes], arrival[nodes])[0])

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
for ax, d, name in [(ax1, d_eff, "effective distance"), (ax2, d_hop, "hop distance")]:
    x = np.array([d[n] for n in nodes], float)
    if ax is ax2:
        x = x + rng.uniform(-0.12, 0.12, x.size)      # 横向抖动，便于看清重叠的点
    ax.scatter(x, arrival[nodes], s=20, alpha=0.6)
    ax.set_xlabel(name); ax.set_title("r = %.2f" % pearsonr([d[n] for n in nodes], arrival[nodes])[0])
ax1.set_ylabel("arrival time (days)")
fig.suptitle("Arrival time follows effective distance, not hops")
plt.tight_layout(); plt.show()

Run在浏览器中运行

首次加载 Python 环境约需 20–40 秒；可以修改参数后重新运行。

集合种群模型与有效距离