基于SEIR 的传染病传播模型

朱 柳 郑婉婷 张贵层

(北方工业大学,北京100000)

纵观历史,一旦出现较大型的传染病传播情况,对世界的经济发展和人民生活都会带来了很大影响。因此想要预测和控制传染病的蔓延,我们就必须了解一般传染病毒的传播机制,定量地研究传染病的传播规律,找出其传播的数学模型。

本文结合实际影响因素建立了一个最接近现实情况的传播模型,并通过建立的预测模型对病毒传播未来一段时间发展态势进行预测。

1 传播模型的建立

1.1 传染病传播机制

针对一般传染病的传播情况,把人群大致分为四类:(1)健康者S(即没有被感染的人群,有被感染的几率)；(2)潜伏者E(已经被感染的人群,但是没有出现症状。根据查阅相关资料[1]潜伏者也具有一定的传染性)；(3)感染者I(即已经被感染并且出现一定的症状)；(4)排除者R(即康复人群和死亡人群,不再参与到感染的过程)。通过相关资料[2],显示潜伏期大约为1—14天,本文对感染者和潜伏者的定义为,两天以下的潜伏期人群称为感染者,三天及以上的潜伏期人群称为潜伏者。

设t 时刻SEIR 四类人口数分别记为s(t)、e(t)、i(t)和r(t),某国的人口总数为N。

α1为感染者I 的单位时间传播率,即单位时间内一个感染者会导致α1S(t)/N 个健康人带上病毒。其中带病毒的人群有k1的比例抵抗力较差从而会直接患病成为感染者I,有k2比例因为自身有较强的抵抗力,没有立即患病,但是在之后也有很大患病几率,成为潜伏者E。

α2为潜伏者E 的单位时间传播率,即单位时间内一个潜伏者会导致α2S(t)/N 个健康人带上病毒,成为潜伏者E,因为潜伏者的传染率比感染者的传染率更弱,大部分潜伏者传染别人使别人变成潜伏者E,而不会使得其直接变成感染者I。

γ1为感染者被严格隔离的比例,γ2为潜伏者中被严格隔离的比例。随着政府干涉以及人民的防范意识增强,隔离率随时间成线性递增关系。

β1为单位时间潜伏者E 变为感染者I 的比例,β2为单位时间潜伏者E 变为健康者S 的比例,变成健康者S 后仍然有被感染的风险。

l1为单位时间内感染者I 的治愈率,l2为单位时间内感染者的死亡率。为直观看出其传播机制,本文做出了病毒的传播图。

图1 传染病传播图

1.2 病毒传播模型微分方程

根据传播机制建立关于s(t),e(t),i(t)和r(t)的方程:

1.2.1 四种人群人数总和为N,保持不变:

1.2.2 等式左边表示单位时间内健康者S 的变化,等式右边从左向右第一部分是正,表示单位时间内潜伏者会有β2比例变成健康者。第二部分是负,表达单位时间未被隔离的感染者(比例为1-γ1)会感染健康者使得健康者变为感染者或潜伏者,传播率为α1。第三部分是负,表达单位时间未隔离的潜伏者(比例为1-γ2)会感染健康者使得健康者变为潜伏者,传播率为α2。

1.2.3 等式左边表示单位时间内潜伏者E 的变化,等式右边从左向右第一部分是正,表示单位时间未被隔离的感染者(比例为1-γ1)会感染健康者使得健康者变为潜伏者,比例为k2,传播率为α1。第二部分是正,表达单位时间未被隔离的潜伏者(比例为1-γ2)会感染健康者使得健康者变为潜伏者,传播率为α2。第三部分为负,表示单位时间内潜伏者会有β1比例变成感染者。第四部分为负,表示单位时间内潜伏者会有β2比例变成健康者。

1.2.4 等式左边表示单位时间内感染者I 的变化,等式右边从左向右第一部分是正,表示单位时间未隔离的感染者(比例为1-γ1)会感染健康者使得健康者变为感染者,比例为k1,感染率为α1。第二部分是正,表达单位时间内潜伏者会有β1比例变成感染者。第三部分为负,表示单位时间内有l1比例的感染者会被治愈成为排除者R,第四部分为负,表示单位时间内有l2比例的感染者会死亡成为排除者R。

1.2.5 等式左边表单位时间内排除者R 的变化,等式右边第一部分为正,表示单位时间内有l1比例的感染者会被治愈成为排除者R,第二部分为负,表示单位时间内有l2比例的感染者会死亡成为排除者R

1.2.6 单位时间内感染者会感染健康人群使得比例为k1的健康者变成感染者,比例为k2的健康者变成潜伏者,它们之和为1。

联立(1)-(6)式进一步整理可得以下方程:

设健康者S、潜伏人群E、感染者I 和康复人群R 初始值分别为s1,e1,i1,r1。

1.3 系数优化模型的建立

通过现有的数据,可以拟合出感染者I 以及排除者R 随时间的变化曲线。在模型中有两个重要的参数:单位时间传播率α1,以及单位时间内潜伏者向感染者转化的比率β1。通过现有资料,大致给出一个范围,即为优化模型的约束条件。通过建立的系数优化模型,可以用智能算法拟合出最符合真实数据的曲线,以得到最佳系数。

2 算例分析

2.1 传播模型系数的估计

人口总N:经过网站[3]查询某国2020 人口总数为143,964,709。

2.1.1 确定系数α1α2

α1:为未被隔离的感染者的传播率,是一个重要参数,对其进行优化得出其值,α2为未被隔离的潜伏者的传播率,由相关文献[4]称,4 成左右的潜伏者有感染能力,在此设α2=0.4×α1,α1范围为:[0.01-1]。

2.1.2 确定系数γ1γ2

γ1γ2分别为感染者的平均隔离率和潜伏者的平均隔离率,没有确定统计的值,隔离率随着时间线性递增,假设:隔离率,γ2=γ1=0.2+0.01×t。

2.1.3 确定系数β1β2

β1为单位时间潜伏者变成感染者的比例,β2为单位时间潜伏者变为健康者的比例由查阅相关资料可得,最终潜伏者中有45%的比例变为感染者,那么55%的比例会变成健康者,基于目前的流行病学调查,潜伏期1-14 天,多为3-7 天,范围为:[0.01-0.45]。

2.1.4 确定系数l1l2

根据真实数据用Matlab2018a 可绘出:

图2 单位时间死亡率和治愈率曲线

记3.29 为第一天,深色线条代表单位时间治愈率,浅色线条代表单位时间死亡率,单位时间死亡率基本保持一样,取其平均值代表其值,算得l2单位时间死亡率为0.001。深色线有一定波动,用cftool 工具进行拟合。得到最拟合方程式。

图3 单位时间治愈率拟合曲线图

l1单位时间治愈率=f(t)=p1×t3+p2×t2+p3×t+p4 p1=2.447e-06 P2=-0.000115 P3=0.001307 P4=0.009866

2.1.5 确定系数k1k2

k1为被感染者感染的健康者变为感染者的比例,k2为被感染者感染的健康者变为潜伏者的概率。潜伏期为1～14 天,在此定义潜伏期2 天以下的就视为感染者,三天及以上视为潜伏者。

根据现有文献研究[4],可以算出潜伏期两天以下的占比约为45%,即k1为45%,潜伏期三天及以上的占比约为k2为55%。

2.2 最佳系数优化解

用MATLAB2018b 编程求解 由粒子群算[5]法拟合出的数据为

α1=0.2923 β1=0.40

拟合值和真实值效果图如图4:

图4 拟合值真实值比较图

2.3 R0 的最终确定

基本传染数(basic reproduction number,R0)是一个流行病学术语,在没有外力介入、所有人都没有免疫力的情况下,一个感染者在具有传染性的这段时间内,平均可以传染给多少人。

在本文的传播模型中,把隔离率设置为0 来模拟没有外力介入的情况。将用粒子群算法拟合出来的单位时间传播率α1=0.2923,和潜伏者转化率β1=0.40 带入到传播模型。

图5 有干涉和无干涉病毒传染比较图

根据相关文献[6]可得:

R0=(1+rTL)(1+rTI)

其中r 表示指数增长的增长率,若用b(t)表示第t 天的新增被感染数,则:

b(t)=b(t-△t)er△t

TL表示感染者潜伏期的平均长度,SI 表示一个感染者被感染的时间和他感染的下一个人被感染的时间的间隔,TI=SI-TL,这三个数据都来自对现实病人情况的观察,通过查询文献[7]可以估计SI 为8.4 天,TL为6 天。

最后剩一个未知数据r 为指数增长率,本文用拟合的无干预(即隔离率为零)数据来拟合r。

所以r=Ln1.1787=0.1644

最终求得R0=(1+rTL)(1+rTI)=(1+0.1644×6)(1+0.1644×2.4)=2.77

3 病毒传播数据的预测

传播模型对病毒传播数据的预测:

用已做好的模型进行预测,记3.29 日为第一天,前面的传播模型有现存感染者I,和R 的数据,用微分方程很容易分别求得累计治愈人数以及累计死亡人数、新增确诊人数为:

当日存在感染者+当日累计治愈+当日累计死亡- (前一日日存在感染者+前一日日累计治愈+前一日累计死亡),用MATLAB 拟合图线(3.29 日为第一天)可得:

图6 病毒传播情况预测图

列出相应表格可得出:

表1 病毒传播情况预测表

同时可以看出在5.14-5.16,将会迎来病毒传播的拐点,感染者人数会下降,新增确诊感染者数也会下降。

4 模型结果分析

通过现有的文献可得一般的R0为2-3,解得的R0大致为2.77,所以传播模型的建立应该具有合理性。通过4.30 日之前的数据,预测出5.1 到5.20 号的数据。根据现有的5.1,5.2 日的真实数据。比较其误差,累出误差比较表可得。(表2)

表2 误差比较表

根据5.1、5.2 日预测和真实值比较可以得出,误差相对小,证明模型合理,结果合理。