人群接触网络中的 SIR 疫情模拟

原创

数据酷客

修改于 2020-05-06 18:05:39

8.5K2

修改于 2020-05-06 18:05:39

文章被收录于专栏：数据科学人工智能数据科学人工智能

查看本案例完整的数据、代码和报告请登录数据酷客（http://cookdata.cn）案例板块。

快速获取案例链接、直播课件：数据酷客公众号内发送“疫情案例”。

视频内容

如何用网络来表示人之间的接触关系？在接触网络中，如何通过 SIR 模型模拟疫情的发展趋势？

本案例将介绍SIR模型，图和网络的基本知识。然后使用 networkx 工具，在生成的随机网络和真实的网络数据上，实现网络中的 SIR 模型进行疫情模拟。

1 SIR 模型介绍

SIR 模型用于计算封闭人群中随着时间推移感染传染性疾病的人数。最早提出来解释在瘟疫（1665-1666年伦敦，1906年孟买）和霍乱（1865年伦敦）等流行病中观察到的感染病人数量的迅速上升和下降。

伦敦大瘟疫是1665年~1666年间，爆发在英国伦敦的大规模传染病，超过8万人死于瘟疫，相当于当时城市人口的1/5。

它假定人口规模是固定的（即无出生和自然死亡等）。传染源的潜伏期是瞬时的，传染持续时间与疾病的长度相同。人群是没有结构的，人之间的接触是完全随机的。恢复之后即获得免疫力，不会继续染病。

1.1 SIR 模型

SIR 是传染病中的最基础最核心的模型，可以使用下面的图形来表示。

如图所示，在 SIR 模型中，人群被划分为以下三类：

易感者（Susceptible）：未感染人群，与感染者接触后可能感染疾病，记为 S(t) 或 S；
感染者（Infectives）：感染人群，感染者可能将疾病传染给易感者，也会以一定概率恢复健康，记为 I(t) 或 I;
恢复者（Recovered）：感染后恢复健康的人，在可能致死的疾病中，也可能包括死亡人群，记为 R(t) 或 R。

SIR 模型包括两个重要的参数：感染率 β，指易感者与感染者接触后染病的概率；恢复率 γ，指感染者恢复的概率。

1.2 使用 Scipy 求解 SIR 模型

上述常微分方程组要直接得到解析解（即把 S,I,R 分别直接写成时间 t 的函数）是相当困难的。我们需要使用数值计算的方法。 Scipy 是 Python 中一个做数值计算的包，我们使用这个包来求解上述方程组。首先，定义一个函数 SIR 来表示 SIR 模型。

def SIR(y,t,beta,gamma):
    S,I,R = y
    dSdt = - S*(I / (S + I + R))*beta
    dIdt = beta*S*I/(S + I + R) - gamma*I
    dRdt = gamma*I
    return [dSdt,dIdt,dRdt]

接下来我们就开始运用 scipy.integrate.odeint() 函数，获得微分方程组的解。

seed = 123 #本案例中将会使用的随机数种子
days = 100 #设置模拟的天数
beta = 0.30 #感染率
gamma = 0.10 #恢复率
N = 150 #人群大小
I0 = 1 #初始感染人数
R0 = 0 #初始恢复人数
S0 = N - I0 - R0

# 设置初始值
y0 = [S0, I0, R0]
from scipy.integrate import odeint
# 求解
solution = odeint(SIR, y0, range(0,days), args = (beta, gamma))

将模拟的结果使用matplotlib工具绘制出来，这里我们直接使用DataFrame对象封装的plot方法。

import matplotlib.pyplot as plt
import pandas as pd
#plt.rcParams['figure.dpi'] = 100
solution_df = pd.DataFrame(solution, columns = ["S","I","R"])
#设置不同人群的显示颜色，易感者为橘色，感染者为红色，恢复者为绿色
color_dict = {"S":"orange","I":"red","R":"green"} 
solution_df.plot(figsize=(9,6),color=[color_dict.get(x) for x in solution_df.columns])

为了更好的演示效果，我们使用pyecharts将图动态地演示出来。

from pyecharts.charts import Line, Grid
import pyecharts
import pyecharts.options as opts
SIR_line = Line().add_xaxis(xaxis_data = solution_df.index)
for col in solution_df.columns:
    SIR_line.add_yaxis(series_name = col,  
             y_axis = solution_df[col].values,  # 添加上证指数数据
             symbol_size = 3,
             color = color_dict[col],
             label_opts = opts.LabelOpts(is_show = False),
             linestyle_opts = opts.LineStyleOpts(width = 1.5), # 设置线宽
             is_smooth = True)
# 在notebook中进行渲染     
SIR_line.render_notebook()

2 网络中的SIR模型

网络可以表示成图 G=(V,E)，其中 V 是节点集合， E 表示边集合，其中每条边由 V 中两个点相连接所构成。如果将人之间的接触关系表示成图，那么图中的节点表示人，边则表示人之间的接触关系。不难想象，如果一个人与他人的接触越多，则在图中该节点与其他节点连接的边也会越多。

在 SIR 模型中，假设人之间是随机接触的。如果人之间的接触关系不是随机的，而是形成了一个接触网络。那么在这个网络中，每个人接触到感染者的概率不再相等，而与他在网络中的位置相关。

与传统 SIR 模型类似，有两个重要的参数：感染率 β 和恢复率 γ。我们需要给每个节点引入一个状态，取值为 S,I,R 中的一种。每一个时间步中，需要动态对每一个节点的状态进行更新。更新规则如下：

如果当前节点是恢复者，则下一步，节点状态依然是恢复者。

如果当前节点是感染者，则下一步，γ 的概率转化为恢复者。1 - γ 的概率依然是感染者。

如果当前节点是易感者，则需要计算其邻居节点中感染者的数量，假设其有 k 个邻居为感染者。则该节点下一步转化为感染者的概率为 1 - (1 - β)k，否则继续保持易感者状态。

对于本次新冠病毒肺炎疫情，如果要使用以上网络的方法对疫情进行模拟，我们还需要构建一个人之间的接触网络。从新闻的疫情报道以及本系列案例中第一次用爬虫抓取的数据，是完全不足以构建接触网络的。需要更多的数据，而这在当前阶段是十分困难的。

本案例中我们采用两种办法简单地构建一个网络结构：使用随机图生成算法生成一个无标度网络；使用一个真实的小型人群接触网络数据集。

3 生成无标度网络进行 SIR 疫情模拟

3.1 无标度网络

统计物理学家把服从幂律分布的现象称为无标度现象，即系统中个体的尺度相差悬殊，缺乏一个优选的标度。于是，满足幂律分布的网络也被称为无标度网络（scale-free network）。

无标度网络中，节点的度 d 满足以下分布：

其中 α 为幂律指数，取值一般在2到3之间。

匈牙利科学家 Barabási 和 Albert 提出了一个 BA 模型（Barabási–Albert model）来生成无标度网络。 BA 模型整体流程如下：

3.2 使用 Networkx 生成无标度网络

Python 中的 Networkx 包提供了方便的随机网络生成函数。其中 barabasi_albert_graph 函数实现了 BA 模型生成无标度网络。主要的参数有网络节点数 m 和新加节点的边数 n 。在我们的场景中，第二个参数的含义是一个人平均与多少人接触。Networkx 包还提供了一系列将网络可视化的函数，能够方便地观察网络的结构。我们这里使用 draw_networkx 函数，它的常用参数包括网络对象、是否显示节点标签（with_labels）、网络的布局（pos）等。

#生成100个节点的BA无标度网络
import warnings
warnings.filterwarnings("ignore")

import networkx as nx #导入networkx包，命名为nx
random_network = nx.barabasi_albert_graph(100,2) # 生成无标度网络，节点数和每个节点边数分别为100和2
#网络可视化
nx.draw_networkx(random_network,with_labels = True,pos = nx.spring_layout(random_network,seed = 1))

3.3 SIR 疫情模拟与分析

节点状态模拟更新函数 updateNodeState

我们使用一个简单的函数来实现一个节点的状态的更新。首先，如果一个节点是恢复者，那么下一步还是恢复者，其节点状态保持不变。如果一个节点是感染者，那么其恢复的概率是 γ。用程序实现的方法为，先均匀生成一个0到1的随机数 p，如果 p < γ，则节点恢复，否则节点依然处于感染状态。

如一个节点是易感者，先要去其邻居节点中看看一共有多少个邻居是感染者，有 k 个邻居是感染者，那么当前节点被感染的概率是 1 - (1 - β)k。我们生成一个0到1的随机数 p，如果 p < 1 - (1 - β)k，则节点被感染，否则不被感染。

updateNodeState 函数实现如下所示，其中输入的参数为网络 G，节点 node，以及 SIR 模型的参数 β 和 γ。

import random
# 根据 SIR 模型，更新节点的状态
def updateNodeState(G,node, beta, gamma):
    if G.nodes[node]["state"] == "I": #感染者
        p = random.random() # 生成一个0到1的随机数
        if p < gamma:   # gamma的概率恢复
            G.nodes[node]["state"] = "R" #将节点状态设置成“R”
    elif G.nodes[node]["state"] == "S": #易感者
        p = random.random() # 生成一个0到1的随机数
        k = 0  # 计算邻居中的感染者数量
        for neibor in G.adj[node]: # 查看所有邻居状态，遍历邻居用 G.adj[node]
            if G.nodes[neibor]["state"] == "I": #如果这个邻居是感染者，则k加1
                k = k + 1
        if p < 1 - (1 - beta)**k:  # 易感者被感染
            G.nodes[node]["state"] = "I"

网络中所有节点状态模拟更新函数 updateNetworkState

有了单个节点状态更新的函数，为了便于后续使用，我们再实现一个整个网络状态进行更新的函数。

def updateNetworkState(G, beta, gamma):
    for node in G: #遍历图中节点，每一个节点状态进行更新
        updateNodeState(G,node, beta, gamma)

为了动态地追踪易感者、感染者和恢复者数量，我们实现一个函数 countSIR 统计三类人群的数量。

# 计算三类人群的数量
def countSIR(G):
    S = 0;I = 0
    for node in G:
        if G.nodes[node]["state"] == "S":
            S = S + 1
        elif G.nodes[node]["state"] == "I":
            I = I + 1
    return S,I, len(G.nodes) - S - I

为了更清楚地演示在网络中疾病的传播过程，我们分别将图中的节点使用不同的颜色进行展示。易感者为橘色，感染者为红色，恢复者为绿色。

def get_node_color(G): #返回每一个节点的颜色组成的列表
    color_list = []
    for node in G:
        #使用我们前面创建的状态到颜色的映射字典 color_dict 
        color_list.append(color_dict[G.nodes[node]["state"]])
    return color_list

在正式开始模拟之前，我们需要进行一些初始化工作，包括节点状态的初始化和 SIR 模型的参数 β 和 γ 的初始化。

ba = nx.barabasi_albert_graph(N,2,seed=seed) 
#初始化节点 state 属性
for node in ba:
    ba.nodes[node]["state"] = "S"
#随机选取一个节点为初始感染者  
ba.nodes[55]["state"] = "I"

在图中开始 SIR 模型的模拟。设置模拟天数，开始执行模拟过程。

# 模拟天数为days，更新节点状态
import matplotlib.pyplot as plt
#fig,ax = plt.subplots(111)
%matplotlib inline
import time
SIR_list = []
for t in range(0,days):
    updateNetworkState(ba,beta,gamma) #对网络状态进行模拟更新
    SIR_list.append(list(countSIR(ba))) #计算更新后三种节点的数量

对模拟的结果进行可视化，查看易感者、感染者和恢复者人数的变化趋势。

df = pd.DataFrame(SIR_list,columns=["S","I","R"])
df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])

我们再次画出在相同的 β 和 γ 参数，相同的 S0,I0,R0 下，传统 SIR 模型的结果。

solution_df.plot(figsize=(9,6),color=[color_dict.get(x) for x in solution_df.columns])

对比两张图，可以看到，在网络中，疫情发展更为迅猛，最高峰的感染人数也远远高于传统 SIR 模型的模拟结果。为什么呢？作为一个开放性的问题，留给大家自己去想吧。

上面的疫情模拟中展示了每一天不同人群的变化，那么在网络中每一天到底是哪些人感染了？我们可以通过 networkx 提供的网络可视化工具深入地分析。通过 networkx.draw_networkx 函数可以方便地将图画出来。

nx.draw_networkx(ba,with_labels = True, node_color = get_node_color(ba), pos = nx.spring_layout(ba,seed = 1))

默认的展示并不是很清晰，我们需要一些额外的设置，例如标签颜色，去掉图片边框，去掉坐标轴刻度，重新设置图的大小等。

fig, ax = plt.subplots(figsize=(9, 6)) #将图的大小设置为 9X6
pos = nx.spring_layout(ba,seed = 1) #设置网络布局，将 seed 固定为 1
ax.axis("off") #关闭坐标轴
plt.box(False) #不显示框
nx.draw(ba,with_labels = True,font_color="white",node_color = get_node_color(ba), edge_color = "#D8D8D8",pos = pos, ax=ax)

将每一天的节点状态以动画演示。我们借助 matplotlib 包中的动画模块 animation 来实现这一效果。

#运行出结果大概需要90秒
#初始化节点 state 属性
for node in ba:
    ba.nodes[node]["state"] = "S"
ba.nodes[55]["state"] = "I"

#绘制网络节点颜色
fig, ax = plt.subplots(figsize=(9, 6))
  
def graph_draw(i,G,pos,ax,beta,gamma): ## 实现动画中每一帧的绘制函数，i为第几帧
    ax.axis("off")
    ax.set_title("day " + str(i) + " 黄色（易感者），红色（感染者），绿色（恢复者）")
    plt.box(False)
    if i == 0: #第一帧，直接绘制网络
        nx.draw(G,with_labels = True, font_color="white", node_color = get_node_color(G), edge_color = "#D8D8D8",pos = pos, ax=ax)
    else: # 其余帧，先更新网络状态，再绘制网络
        updateNetworkState(G, beta, gamma)
        nx.draw_networkx_nodes(G, with_labels = True, font_color="white", node_color = get_node_color(G),pos = pos, ax=ax) 
    plt.close()
    
#演示网络动态变化
import matplotlib.animation as animation
from IPython.display import HTML

animator = animation.FuncAnimation(fig, graph_draw, frames= range(0,days),fargs=(ba,pos,ax,beta,gamma),  interval=200)
HTML(animator.to_jshtml())

3.4 初始感染者的影响

初始感染者在网络中位置不一样，对疾病的传播也会造成影响。这里我们用两种策略进行对比分析：

随机选择初始感染者；
选择度数最高的点作为初始感染者。

首先，我们随机选择一个种子节点。

import numpy as np
seed_node = np.random.randint(0,N)
print(seed_node)
#初始化节点 state 属性
for node in ba:
    ba.nodes[node]["state"] = "S"
ba.nodes[seed_node]["state"] = "I"

# 模拟days天，更新节点状态
SIR_list = []
for t in range(0,days):
    updateNetworkState(ba,beta, gamma)
    SIR_list.append(list(countSIR(ba))) 
df = pd.DataFrame(SIR_list,columns=["S","I","R"])
df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])

然后，计算图中节点的度，选择度最高的节点作为种子感染者。

# 使用networkx 的 degree_centrality 函数计算图中节点的度中心度
node_degree = nx.degree_centrality(ba)

node_degree_df = pd.DataFrame.from_dict(node_degree, orient = "index", columns=["degree"])
node_degree_df = node_degree_df.reset_index().rename(columns = {"index":"node"})

#查看度数最高的节点
node_degree_df.sort_values(by = "degree",inplace = True,ascending= False)
node_degree_df.head()

下面我们将度数最高的节点设置为种子感染者，观察疫情变化曲线。

seed_node = int(node_degree_df.values[0,0])
print(seed_node)
#初始化节点 state 属性
for node in ba:
    ba.nodes[node]["state"] = "S"
ba.nodes[seed_node]["state"] = "I"

# 模拟days天，更新节点状态
SIR_list = []
for t in range(0,days):
    updateNetworkState(ba,beta, gamma)
    SIR_list.append(list(countSIR(ba))) 
df = pd.DataFrame(SIR_list,columns=["S","I","R"])
df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])

4 真实人群网络中的 SIR 疫情模拟

我们这里选用一个真实的人与人之间的接触网络来模拟疫情的传播。Infectious 是 2009 年在都柏林举办的主题为 “INFECTIOUS: STAY AWAY” 的展览活动中人与人之间面对面接触的网络。图中一共包含 410 个节点，每个节点表示一个人。如果两个人之间有超过20秒以上的面对面接触，则它们之间存在一条边。原始数据集中两个节点之间可能存在多条边，为了简化分析我们只保留其中的一条边。数据集来源于网站KNOECT。原始图数据见文件out.sociopatterns-infectious.csv，去除了多余边的网络文件见infectious.csv。

首先，我们使用 networkx 的 read_edgelist 函数加载网络，并简单可视化。

infectious_network = nx.read_edgelist("./input/infectious.csv",delimiter=",")
fig, ax = plt.subplots(figsize=(24, 16)) #节点较多，将图片大小也调整大些
pos_infectious = nx.spring_layout(infectious_network,seed = 22)
ax.axis("off")
plt.box(False)
nx.draw(infectious_network,with_labels = True,font_color="white", node_color = "orange", edge_color = "#D8D8D8",pos = pos_infectious, ax=ax)

选择最上图中左边的 397 号节点作为种子节点，观看疫情的变化。注意，398 号节点处于网络中边缘的位置,在模拟过程中有可能尚未将疾病传播出去 397 号就恢复了健康，因此疾病不会在网络中继续传播。

beta = 0.10 # 为了更好观察，我们减小传染率参数
gamma = 0.05
#初始化节点 state 属性
for node in infectious_network:
    infectious_network.nodes[node]["state"] = "S"
infectious_network.nodes["397"]["state"] = "I"

# 模拟days天，更新节点状态
SIR_list = []
for t in range(0,days):
    updateNetworkState(infectious_network, beta, gamma)
    SIR_list.append(list(countSIR(infectious_network)))
    
df = pd.DataFrame(SIR_list,columns=["S","I","R"])

df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])

选取 272 号节点作为种子节点，观察疫情发展。

#初始化节点 state 属性
for node in infectious_network:
    infectious_network.nodes[node]["state"] = "S"
infectious_network.nodes["272"]["state"] = "I"

# 模拟days天，更新节点状态
SIR_list = []
for t in range(0,days):
    updateNetworkState(infectious_network, beta, gamma)
    SIR_list.append(list(countSIR(infectious_network)))
    
df = pd.DataFrame(SIR_list,columns=["S","I","R"])
df.plot(figsize=(9,6),color=[color_dict.get(x) for x in df.columns])

观察上述两个图中疫情高峰到来的时间，可以看到网络中重要的节点可以更快地传播疾病。这些人也是疫情防控中的关键。最后我们通过一个动画来感受下在这个真实的接触网络中疫情的传播态势。

#初始化节点 state 属性
#plt.rcParams['figure.dpi'] = 150
plt.rcParams['animation.embed_limit'] = 217277560

for node in infectious_network:
    infectious_network.nodes[node]["state"] = "S"
infectious_network.nodes["397"]["state"] = "I"

fig, ax = plt.subplots(figsize=(16, 12))
animator = animation.FuncAnimation(fig, graph_draw, frames= range(0,days),fargs=(infectious_network,pos_infectious,ax,beta,gamma),  interval=200)
HTML(animator.to_jshtml())