SARS疾病传播的数学模型.docx
《SARS疾病传播的数学模型.docx》由会员分享，可在线阅读，更多相关《SARS疾病传播的数学模型.docx（20页珍藏版）》请在冰点文库上搜索。
SARS疾病传播的数学模型
SARS传播的数学模型
摘要
本文针对SARS的传播建立了数学模型。
首先，对附件1提供的早期模型，认为“传染概率”的说法欠妥，传染期限L的确定缺乏医学上的支持，使模型的说服力降低。
模型中借鉴广东香港的参数来预测北京的疫情走势，不失为一种方法，但在不同地区因政策，地域的不同，病毒的传播和控制呈现不同的特点，使不同城市之间的可比性降低。
故借鉴法存在一定的适用范围，且不能对首发城市进行预测。
对于第二问，在分析常用传染病模型的局限性后，文中把患者所处的状态明确划分为潜伏阶段、发病阶段和隔离阶段，根据各阶段的转化关系建立了第一个数学模型。
考虑到发病和被隔离等事件发生的随机性，本文在原有模型的基础上适当改进，建立了随机模拟模型。
通过对5月10日以前数据的拟合，并经过500次模拟，对北京的疫情进行了预测：
7月上旬北京将基本解除疫情，累计病例约2800多人。
预测结果与实际情况符合得很好。
另外，改变有关参数，发现提前5天采取严格的隔离措施，将使疫情解除的时间提前约10天，累计人数降至1958人；若延迟5天采取措施，疫情将推迟11天，累计人数达4487人。
根据这些预测，文中对卫生部门采取控制措施提出了相关建议。
对第三个问题，本文研究SARS对入境旅游人数的影响，建立了数学模型。
通过数据拟合的方法确定日增长病例数对旅游人数的影响，预测9~12月份入境旅游人数分别为24.02，36.06，33.04，25.85万人。
与往年同期相比，9月降低了23.5个百分点，10月以后影响逐步减小，经济进入恢复时期。
对于第四个问题，给报刊写了一篇通俗短文，说明了建立传染病数学模型的重要性。
最后在模型的评价中，对该模型优于原附件1模型的方面作了说明，特别说明了建立一个真正能预测和为预防、控制提供可靠、足够的信息的模型需要满足的条件和困难之处。
一、问题的提出
2002年至2003年，SARS（严重急性呼吸道综合症,俗称非典型肺炎）悄然无息地靠近我们的生活，在潜伏一段时间后忽然爆发，在全球掀起了轩然大波。
作为重灾区的国家之一，我国的经济发展和人民生活受到了很大的影响。
我们从中得到了许多重要的经验和教训，认识到定量研究传染病的传播规律、为预测和控制传染病蔓延创造条件的重要性。
对此，要求对SARS的传播建立数学模型，具体要求如下：
1、对附件1所提供的一个早期的模型，评价其合理性和实用性。
2、对SARS的传播建立一个自己的模型，并说明：
（1）为什么优于附件1中的模型；
（2）怎样才能建立一个真正能够预测以及能为预防和控制提供可靠、足够的信息的模型，以及这样做的困难之处。
（3）对于卫生部门所采取的措施做出评论，如：
提前或延后5天采取严格的隔离措施，对疫情传播所造成的影响做出估计。
（附件2提供的数据供参考。
）
3、收集SARS对经济某个方面影响的数据，建立相应的数学模型并进行预测。
（附件3提供的数据供参考。
）
4、给当地报刊写一篇通俗短文，说明建立传染病数学模型的重要性。
二、对早期模型的评价
附件1的模型主要采用“数据拟合”和“借鉴参数”的方法对北京疫情走势进行预测。
在数据拟合方面，该模型中有两个疑点：
1、感染期限L的确定。
由于被严格隔离、治愈、死亡等原因，感染者在某一时段后不再具有对易感人群的传染力，故对病毒的传染加上感染期限是合理的。
但在对该参数的确定上，作者为了较好地拟合各阶段的数据,通过人为调试来确定L的取值，缺乏医学上的支持，使模型的说服力减弱，合理性和可靠性大大降低。
2、文中认为“K代表某种环境下一个人传染他人的平均概率”。
但从模型的公式中可以看出，参数K的实际意义是一个病人平均每天传染其他人的个数。
两者之间有实质的区别，文中的说法显然不妥。
从预测思想来看，该模型是借鉴先发地区——广东、香港的有关参数对北京的疫情进行预测的。
由于广东、香港的疫情和控制都在北京之前，已经过了高峰期，到5月8日为止每日新增病例已降至10来例，基本处于后期控制阶段。
而当时北京的疫情刚过了高峰期，正处于社会剧烈调整时期，数据较为凌乱，略有下降趋势，但不明显。
可见在当时，采取这种借鉴是无奈之举。
但是由于城市之间的政策，风俗习惯等不同，城市之间的可比性不强，借鉴存在很大的局限性。
如在香港，由于对传播机制认识不足，中途又出现高度感染的特殊情况。
另外使用借鉴法无法对首发城市进行预测。
三、传播模型
（一）问题的分析
在SARS爆发的初期，由于潜伏期的存在，人们对病毒传播速度和危害程度的认识不够，未能及时识别这一传染病的存在。
但当病患数不断增加，政府开始采取应对措施对其进行控制，同时社会舆论加大宣传力度，人们的警觉性提高，病毒的传播速度下降。
因此，我们通常把传染病的传播模式近似分为两个阶段：
第一、自由传播阶段（即控前阶段）：
在采取切实有效的控制措施之前的一段时间。
第二、控后阶段：
介入人为因素之后的一段时间。
由于SARS的传播涉及的因素很多，如潜伏期、人群的迁入迁出，感病者的数量、易感者的数量、传染率和治愈率的大小等。
而且在以上因素中，潜伏期的大小、传染率和治愈率的大小因人而易，具有一定的随机性。
不可能一开始就把所有的因素全部考虑在内建立模型。
对此，我们将作出相应的假设进行简化。
分析附件2所给出的数据，发现6月1日至15日，已确症病例累计数为2522人，其中夹杂3天累计数为2523人。
但6月16日后累计数降至2521人，认为累计数的减少可能是误诊引起的。
由于误诊的可能性很低，故在这里忽略不计。
（二）基本假设
1、国家卫生部提供的北京疫情统计真实可信。
（误诊数仅为1，可忽略不计）。
2、由于非典的主要传播途径是近距离接触，通过受感染者咳嗽或打喷嚏时产生的飞沫传播，这里将所有传播途径都视为与病源的直接接触。
3、不考虑出生与自然死亡的过程和人群的迁入迁出（或认为迁入和迁出基本平衡），认为疾病传播期间所考察地区的总人数为常数。
4、根据国家卫生部资料可知处于潜伏期的SARS病人不具有传染性。
5、目前尚不清楚康复患者是否具有免疫力，但据国家卫生部资料可知康复后的病患无一例复发。
故假设康复患者退出传染系统。
6、根据资料显示，SARS病毒的潜伏期一般为2~7天，平均约为5天。
（这一条件将在后期的模型中有所改动）
（三）常用基本模型
目前常用的传染病模型，通常将传染病流行范围内的人群分为三类：
S类：
易感者，指未得病者，但与感病者接触后容易受到感染。
I类：
感病者，指染上传染病的人。
R类：
移出者，指因患病而被隔离，或因病愈而具有免疫力的人，他们即非感病者，也非易感者，实际上他们已经退出了传染病系统。
并通过三类之间的互相转化关系建立微分方程组进行求解：
（1）
变量和符号说明：
——传染率：
每个病人平均每天有效接触（足以使被解除者感染）的人数。
——退出率：
单位时间内治愈和死亡人数占感病者人数的百分数。
S（t）——易感人群的总数。
I（t）——感病者总数。
R（t）——退出者总数。
N——一个城市总人口数。
观察附件二中给出的数据，我们发现截至6月23日，感病者累计为2521人，远远小于北京城市的总人口数150万人，故认为感病者和退出者对易感人群的总数影响不大，易感者总人数I为一常数。
原方程变形为：
（2）
注意到退出者不是我们研究的范围，故方程组
（2）实际上是一个常微分方程
（3）
其中
，
不难用分离变量法解出：
（4）
其中I0为初始值。
根据以上分析我们可以看出，常微分方程的传染病模型只适用于病例数与总人口数具有可比性的情况。
当病例数远小于总人口数时，常微分方程模型的实质与附件1的模型相同，感病人数将随时间以指数增长。
考虑这一特点，我们用计算机跟踪病毒的个体传播情况，建立了模拟模型。
（四）计算机模拟模型：
在该模型中，我们将传染系统中的人分为五类：
自由携带者（
）——身上携带病毒并均匀散布在人群中的患者，根据基本假设自由携带者在潜伏期内不具有传染力，
日增患者（
）——每天被医疗部门发现并加以隔离的感病者
被隔离者（
）——因曾与自由携带者接触而被怀疑携带SARS病毒的人
有效接触者（
）——每日与自由携带者接触并感染上病毒的人
无效接触者（
）——每日与自由携带者接触但未染上病毒的人
并作出如下假设：
1、由于传染性SARS最初（1~2天）的症状通常为发热（
），发热通常为高热[1]。
症状明显，易于辨认，故可认为自由携带者发病后当天或第二天就立即入院治疗，入院后不会再参与疾病的传播。
2、根据实际情况，假设SARS病人被发现的三天内，有关部门将采取措施，将部分与病源有效接触者隔离，这部分人即使发病后也不会参与疾病的传播。
3、与病源有效接触者必然发病。
根据基本假设，潜伏期一般为2至7天，这里取为5天。
（这一假设在改进模型中有进一步的讨论。
）
另外，对模拟模型中出现的符号变量说明如下：
——有效接触率，表示一个自由携带者平均每天有效接触的人数。
——无效接触率，表示一个自由携带者平均每天无效接触的人数。
——与自由携带者接触后（包括有效接触和无效接触）的人群中可以控制的人数所占的百分比。
模拟模型中个体传播情况如图1所示：
自由携带者1~5天处于潜伏期，不具有传染能力；5天后发病，发病后每天有效接触
人，2天后（第7天）被隔离，在隔离前每天无效接触
人。
与病源接触后的可控人群（占接触者总人数的
）在3天后被视为疑似病人。
疑似病人中的有效接触者在接触病源的第7天被发现确认为日增患者，而有效接触者中其他人作为自由携带者留在人群中，继续这之前的个体传播。
图1：
个体传播示意图
用数学模型描述各个变量之间的关系如下：
（5）
由于北京在4月20日才开始建立每日疫情报道制度，故认为政府采取严格的隔离措施开始于4月20日，以这一天为分界线，之前属于自由传播阶段。
根据这一模型，用计算机模拟北京5月10日之前SARS的传播情况，并对5月10日以后的传播情况进行预测。
图2：
5月10日以前数据拟合图
图3：
5月10日以后的预测曲线
通过图2两条曲线的拟合，得到控前的有效接触率（表征病毒的传染力）k=1.35>1,可控率（表征政府的控制力度）
；控后k1=0.8，
=0.7。
根据这两个参数作出5月10日的预测曲线见图3。
，根据预测，北京将在第97天（6月下旬）实现零增长，累计病例数2448人。
经过分析，以上确定型的模拟模型存在以下两点问题：
第一、SARS病毒的潜伏期一般是2~7天，模型将潜伏期确定为5天，从感染到被发现的时间确定为7天，可以明显地看到以7天为传播周期的曲线变动，而在实际曲线中，虽然数据有上下波动的趋势，但周期性并不明显。
第二、在采用隔离政策时，模型假定与病源接触的人群以固定的比例
受到控制，这种假设加大了人为主观因素的影响。
基于以上两点，在拟合5月10日之前的数据时，得到的有效接触率与实际统计数据有所偏差。
这种偏差降低了模型的可信度。
基于此，我们查找有关统计数据，为参数的确定寻求医学上的支持，并以随机模拟取代完全确定性的模拟，对原模型进行改进，建立了随机模拟模型。
（五）改进后的随机模拟模型
1、根据上文的分析，在原有模型的基础上作出如下改进假设：
假设一：
假设潜伏期的长短服从以5天为期望值，2天为方差的正态分布，即认为自由携带者将以68%的概率在3~7天内发病，以95%的概率在1~9天内发病。
假设二：
假设与病源接触的人群每天均以一定的比例
被隔离。
该比例服从以
为期望值，0.05为方差的正态分布。
2、参数的确定：
控前k1——5月23日，美国《科学》杂志及其网站上发表了两篇有关SARS的流行病学论文及一篇长篇评论[2]，研究确定SARS病患的基本传播率（即文中的有效接触率）为3。
控前
——表示有效接触者中在发病前未能受到有效隔离的人所占的比例。
由于在SARS传播前期，对疾病未能作良好的监控，而在潜伏期内病患无法得知自己的真实情况。
所以在这一阶段不妨设
，即几乎没有病患在潜伏期内被隔离。
这与控前的实际情况是吻合的。
控后k1——控后的有效接触率实际上受多种因素综合影响。
比如人们不愿外出，减少了与病源的接触机会，更注意卫生等，这些因素均使有效接触率明显降低。
英国和香港的一个研究小组的研究表明，传染率（即本文的有效接触率）控后仍有较长一段时间在1左右徘徊。
为简化模型，我们取k1=1。
控后
——控后的隔离措施被细化为几个等级，措施的级别越高，控制力度越大，社会付出代价越高，难度越大[3]。
经研究发现，采取隔离措施后短期内使未隔离的病患在总自由携带者中占40%~90%的比例[3]，即
在0.4在0.9之间。
由于各地的控制力度不同，对应的
不同，对北京4月20日以前的数据进行拟合后得到当
时，拟合较好。
从中可以看出北京的措施是比较严厉的。
3、预测结果
由于缺乏4月20日之前的数据，不妨设初始6天的病例数分别为1，2，3，4，5，6（之后将对初始值的合理性进行讨论）。
为了较好地拟合数据，采用3月15日作为病毒传播的开始时间，认为在病毒传染前期，北京主要以输入型传播为主，传播机理不明，输入病例数量无法确定。
采用以上数据和相关参数，通过计算机编程，产生正态分布的随机数，并对传染情况进行500次模拟，取其平均数得到预测结果见表1。
表1：
实际值与预测值的比较
日期
5.11
5.12
5.13
5.14
5.15
5.16
5.17
5.18
实际值
38
39
43
23
18
17
15
14
预测值
43.
40
37
30
30
26
23
20
日期
5.19
5.20
5.21
5.22
5.23
5.24
5.25
5.26
实际值
3
7
0
12
9
25
9
5
预测值
18
16
15
14
13
11
10
9
日期
5.27
5.28
5.29
5.30
5.31
6.1
6.2
6.3
实际值
8
2
3
3
1
1
0
0
预测值
9
8
7
6
6
5
5
4
日期
6.4
6.5
6.6
6.7
6.8
6.9
6.10
6.11
实际值
0
0
0
0
0
0
0
0
预测值
4
3
3
3
2
2
2
2
图4：
改进后的拟合曲线
图5：
改进后的预测曲线
图6：
提前5天采取措施的预测曲线
若卫生部门提前5天采取严格的隔离措施，预测将会出现图4的情况：
累计病例总数为1958人，到第100天就可实现零增长。
若滞后5天才采取严格的隔离措施，预测曲线如图5所示，疫情将延迟至第121天才能基本消除，届时累计病例总数将达到4487人。
若卫生部门的隔离措施滞后15天才得到全面实施，将出现图6的情况，到第139天疫情才能基本消除，累计病例数将增至10208人,比原来增长了约8000多名患者。
图7：
滞后5天的预测曲线
图8：
滞后15天的预测曲线
若社会及政府等相关部门加大宣传力度，人们警觉性提高，使有效接触率降低，设k=0.8，得到预测曲线图7，预计到第92天疫情解除，累计病例2395人。
若卫生部门实施严格隔离措施，加大对与自由携带者相接触的人群的控制力，设可控率提高到
，
。
得到预测曲线图8，预计到第80天疫情解除，累计病例2545人。
表2：
改变有关参数后的预测值
改变相关参数
累计病例数
实现零增
长的时间
与预测病例
数的差值
与预测实现零增长的日期相比
最初预测
2816
109
提前5天
1958
100
-858
-9
滞后5天
4487
121
1671
+12
滞后15天
10208
139
7392
+30
加大控制力度
2545
80
-271
-29
加大宣传力度
2395
92
-421
-18
以上预测结果表明，加大宣传力度，可以明显减少感染总数；而加大控制力度，可以显著减少疫情延续的时间。
根据以上预测结果，对SARS疫情的控制提出以下建议：
1、应尽早采取有效的控制措施。
提早采取严格的隔离措施，将会使疫情提早得到控制或解除。
2、隔离和宣传并重。
一方面加大控制力度，提高接触病源的人群的可控率，提前结束疫情；另一方面加大宣传力度，提高全民健康意识，切断传染途径，从而减少自由携带者的有效接触人数。
这样即使防范措施晚了一些，出现了高峰期，仍然能够在一定的时间内解除疫情。
四、模型的评价和改进方向
（一）模型的评价：
1、附件1的模型中，未能明确指出自由携带者和已确症患者之间的区别，未能考虑到潜伏期不具备传染力的情况。
而在本文中的两个模拟模型，对各个时期的带菌者作了严格的界定，认为只有未受控（未被隔离及入院就诊）的带菌者在病发后的1~2天内才具有传染能力，之后会立即入院治疗。
符合实际情况，揭示了传染的主要途径和传播机理。
2、改进后的模拟模型对潜伏期和带菌者被发现的时间按照实际规律，作了正态模拟，具有明确的实际意义。
3、模拟模型是完全依照北京市5月10日以前的数据拟合得到的，不需要借鉴其他城市的数据，可以预测首发城市的疫情发展，且结果与实际情况符合的很好。
4、正是由于我们模型逻辑的准确性，如果经过深入调查获取实际数据，再对模型进行一定的修改，就可立即用于建立前期的预测模型，而不需要等待疾病的发展来获得大量数据进行拟和。
（二）理想模型建立的方案及困难
要真正建立能够预测以及能为预防和控制提供可靠，足够的信息的模型。
要怎样做，以及困难。
事实上，真正情况要比我们的模型复杂的多，比如，由于初期传播机理不明，医护人员为收到好的保护，很多受到传染。
而传播率事实上也是按照某种趋势变化的，他与人们的警惕心理等因素有关。
要建立这样的模型，我们认为：
1、要不断的研究传染病理，完善模型。
2、建立实时反馈机制，不断修改模型参数，以求更准确的预测以后的数据。
3、尽量以实际的研究数据建模，力求准确反映事件的实质。
这样，当人为修改参数进行模拟预测时，可以为采取什么措施效果最好做出较准确的建议。
4、研究各种措施实际对疾病传播的影响。
要做到这些，也恰恰是建模的难点，因为当涉及到新的传染病的时候，以上数据或参数都必须通过长期的研究或拟合才能得到的。
比如措施对疾病传播的影响，是要统计大量实际数据才能获得的。
以上也可以看作我们模型的改进方向。
五、参考文献
[1]范学工龙云铸，严重急性呼吸综合征——非典型性肺炎讨论，中国现代医学杂志，第13卷第8期，1-4页，2003年4月
[2]济南市卫生局/机构设置/疾病控制处，SARS传染能力及超级传染者揭秘，传染性非典型肺炎防治，2003-6-2
[3]中国科学院研究生院管理学院“SARS”管理与控制研究小组，对各种措施的效果分析和结论，研究简报，第1期
附件
确定型模拟程序
functionlogi（N）
%default:
N=40
%脚本程序：
（下一行）
%N=40;logi（N）;t=a（find（a>0.1））,sum（t）,plot（1:
length（t）,t,'g'）
clear
f=zeros（1,200）;%f表示自由的带菌者数
a=zeros（1,200）;%a表示新增病例数
b=zeros（1,200）;%b表示被隔离的人数
%f（1:
7）=[60,70,80,100,143,106,105];
f（1:
7）=[1,2,3,4,5,6,7];
%f
（1）的起点：
beta=.5;k=1.35;T=5;l=.01;
fori=1:
N
a（i+6）=f（i）+a（i+6）;
b（i+8）=（7*T*l+2*k）*beta*f（i）/3+b（i+8）;
b（i+9）=（7*T*l+2*k）*beta*f（i）/3+b（i+9）;
b（i+10）=（7*T*l+2*k）*beta*f（i）/3+b（i+10）;
a（i+11）=2*beta*k*f（i）/3+a（i+11）;
a（i+11）=2*beta*k*f（i）/3+a（i+11）;
a（i+11）=2*beta*k*f（i）/3+a（i+11）;
f（i+5）=（1-beta）*k*f（i）+f（i+5）;
f（i+6）=（1-beta）*k*f（i）+f（i+6）;
end
beta=.7;k=.8;
fori=N:
180
a（i+6）=f（i）+a（i+6）;
b（i+8）=（7*T+2*k）*beta*f（i）/3+b（i+8）;
b（i+9）=（7*T+2*k）*beta*f（i）/3+b（i+9）;
b（i+10）=（7*T+2*k）*beta*f（i）/3+b（i+10）;
a（i+11）=2*beta*k*f（i）/3+a（i+11）;
a（i+11）=2*beta*k*f（i）/3+a（i+11）;
a（i+11）=2*beta*k*f（i）/3+a（i+11）;
f（i+5）=（1-beta）*k*f（i）+f（i+5）;
f（i+6）=（1-beta）*k*f（i）+f（i+6）;
end
%N=40;logi（N）;t=a（find（a>0.1））,sum（t）,plot（1:
length（t）,t,'g'）
AA=[60,70,80,100,143,106,105,81,103,111,126,85,148,93,113,83,105,62,94,63,89,87,41,50,38,...
39,43,23,18,17,15,14,3,7,0,12,9,25,9,5,8,2,3,3,1,1,0,0,0,0,0,1,-1,0,0,1,0,-1,0,0,-1,0,0,...
0,0,0,0,0];
figure;
plot（37:
36+length（AA）,AA,'r:
'）;hold;
plot（1:
100,a（1:
100）,'b--'）;hold
legend（'实际','拟和'）
%以a<0.1结束
%N=40:
2448--------------103天
%提前5天：
1811----------------97天
%滞后5天：
3274----------------111天
%滞后15天：
5741------------125
%差值：
1462
升级会员 