2003CMCM
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来例,基本处于后期控制阶段。而当时北京的疫情刚过了高峰期,正处于社会剧烈调整时期,数据较为凌乱,略有下降趋势,但不明显。可见在当时,采取这种借鉴是无奈之举。
但是由于城市之间的政策,风俗习惯等不同,城市之间的可比性不强,借鉴存在很大的局限性。如在香港,由于对传播机制认识不足,中途又出现高度感染的特殊情况。另外使用借鉴法无法对首发城市进行预测。
1
三、 传播模型
(一)问题的分析
在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类:移出者,指因患病而被隔离,或因病愈而具有免疫力的人,他们即非感病者,也非易感者,实际上他们已经退出了传染病系统。
并通过三类之间的互相转化关系建立微分方程组进行求解:
2
?dS?dt??kIS??dI?kIS?hI??dt?dR?hI??dt??S?I?R?N (1)
变量和符号说明:
k ——传染率:每个病人平均每天有效接触(足以使被解除者感染)的人
数。
h ——退出率:单位时间内治愈和死亡人数占感病者人数的百分数。
S(t)——易感人群的总数。
I(t)——感病者总数。 R(t)——退出者总数。
N——一个城市总人口数。
观察附件二中给出的数据,我们发现截至6月23日,感病者累计为2521人,远远小于北京城市的总人口数150万人,故认为感病者和退出者对易感人群的总数影响不大,易感者总人数I为一常数。原方程变形为:
?dI?kIN?hI??dt (2) ??dR?hI??dt注意到退出者不是我们研究的范围,故方程组(2)实际上是一个常微分方程
dIdt?kNI?hI??I (3)
其中??kN?h,
不难用分离变量法解出:
I(t)?I0e??t (4)
其中I0为初始值。
根据以上分析我们可以看出,常微分方程的传染病模型只适用于病例数与总人口数具有可比性的情况。当病例数远小于总人口数时,常微分方程模型的实质与附件1的模型相同,感病人数将随时间以指数增长。
考虑这一特点,我们用计算机跟踪病毒的个体传播情况,建立了模拟模型。
(四)计算机模拟模型:
在该模型中,我们将传染系统中的人分为五类:
自由携带者(f[t])——身上携带病毒并均匀散布在人群中的患者,根据基
3
本假设自由携带者在潜伏期内不具有传染力,
日增患者(x[t])——每天被医疗部门发现并加以隔离的感病者
被隔离者(y[t])——因曾与自由携带者接触而被怀疑携带SARS病毒的人 有效接触者(z1[t])——每日与自由携带者接触并感染上病毒的人 无效接触者(z2[t])——每日与自由携带者接触但未染上病毒的人 并作出如下假设:
1、由于传染性SARS最初(1~2天)的症状通常为发热(?38o),发热通常为高热[1]。症状明显,易于辨认,故可认为自由携带者发病后当天或第二天就立即入院治疗,入院后不会再参与疾病的传播。
2、根据实际情况,假设SARS病人被发现的三天内,有关部门将采取措施,将部分与病源有效接触者隔离,这部分人即使发病后也不会参与疾病的传播。
3、与病源有效接触者必然发病。根据基本假设,潜伏期一般为2至7天,这里取为5天。(这一假设在改进模型中有进一步的讨论。)
另外,对模拟模型中出现的符号变量说明如下:
k1 ——有效接触率,表示一个自由携带者平均每天有效接触的人数。 k2 ——无效接触率,表示一个自由携带者平均每天无效接触的人数。
(包括有效接触和无效接触)的人群中可以控制?——与自由携带者接触后
的人数所占的百分比。
模拟模型中个体传播情况如图1所示:自由携带者1~5天处于潜伏期,不具有传染能力;5天后发病,发病后每天有效接触k1人,2天后(第7天)被隔离,在隔离前每天无效接触k2人。与病源接触后的可控人群(占接触者总人数的?)在3天后被视为疑似病人。疑似病人中的有效接触者在接触病源的第7天被发现确认为日增患者,而有效接触者中其他人作为自由携带者留在人群中,继续这之前的个体传播。
有效接触者 自由携带者 第6,7天 3天后 自由携带者 无效接触者 被隔离者 7天中 7 天后 日增患者 3天后 4天后 日增患者 被隔离者
4
图1:个体传播示意图
用数学模型描述各个变量之间的关系如下:
?z1[i?5]?k1f[i]?z[i?6]?k1f[i]?1??z2[i?j]?k2f[i](j?0,1,2,?6) (5) ?f[i]?(1??)z[i]1???x[i]?f[i?6]??z1[i?6]由于北京在4月20日才开始建立每日疫情报道制度,故认为政府采取严格
的隔离措施开始于4月20日,以这一天为分界线,之前属于自由传播阶段。
根据这一模型,用计算机模拟北京5月10日之前SARS的传播情况,并对5月10日以后的传播情况进行预测。
图2:5月10日以前数据拟合图
5
图3:5月10日以后的预测曲线
通过图2两条曲线的拟合,得到控前的有效接触率(表征病毒的传染力)k=1.35>1,可控率(表征政府的控制力度)??0.5;控后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天内发病。
假设二:假设与病源接触的人群每天均以一定的比例?被隔离。该比例服从以1??为期望值,0.05为方差的正态分布。
2、参数的确定:
控前k1——5月23日,美国《科学》杂志及其网站上发表了两篇有关SARS
6
的流行病学论文及一篇长篇评论[2],研究确定SARS病患的基本传播率(即文中的有效接触率)为3。
控前?——表示有效接触者中在发病前未能受到有效隔离的人所占的比例。由于在SARS传播前期,对疾病未能作良好的监控,而在潜伏期内病患无法得知自己的真实情况。所以在这一阶段不妨设??0.90,即几乎没有病患在潜伏期内被隔离。这与控前的实际情况是吻合的。
控后k1——控后的有效接触率实际上受多种因素综合影响。比如人们不愿外出,减少了与病源的接触机会,更注意卫生等,这些因素均使有效接触率明显降低。英国和香港的一个研究小组的研究表明,传染率(即本文的有效接触率)控后仍有较长一段时间在1左右徘徊。为简化模型,我们取k1=1。
控后?——控后的隔离措施被细化为几个等级,措施的级别越高,控制力度越大,社会付出代价越高,难度越大[3]。经研究发现,采取隔离措施后短期内使未隔离的病患在总自由携带者中占40%~90%的比例[3],即?在0.4在0.9之间。由于各地的控制力度不同,对应的?不同,对北京4月20日以前的数据进行拟合后得到当??0.4时,拟合较好。从中可以看出北京的措施是比较严厉的。
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
7
图4:改进后的拟合曲线
图5:改进后的预测曲线
8
图6:提前5天采取措施的预测曲线
若卫生部门提前5天采取严格的隔离措施,预测将会出现图4的情况:累计病例总数为1958人,到第100天就可实现零增长。
若滞后5天才采取严格的隔离措施,预测曲线如图5所示,疫情将延迟至第121天才能基本消除,届时累计病例总数将达到4487人。
若卫生部门的隔离措施滞后15天才得到全面实施,将出现图6的情况,到第139天疫情才能基本消除,累计病例数将增至10208人,比原来增长了约8000多名患者。
图7:滞后5天的预测曲线
9
图8:滞后15天的预测曲线
若社会及政府等相关部门加大宣传力度,人们警觉性提高,使有效接触率降低,设k=0.8,得到预测曲线图7,预计到第92天疫情解除,累计病例2395人。
若卫生部门实施严格隔离措施,加大对与自由携带者相接触的人群的控制力,设可控率提高到??0.7,??0.3。得到预测曲线图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、隔离和宣传并重。一方面加大控制力度,提高接触病源的人群的可控率,提前结束疫情;另一方面加大宣传力度,提高全民健康意识,切断传染途径,从而减少自由携带者的有效接触人数。这样即使防范措施晚了一些,出现了高峰期,仍然能够在一定的时间内解除疫情。
10
四、 SARS影响经济模型
(一)问题的分析
从题中数据可以看出SARS对旅游行业造成了严重的影响,影响最大的是五月,海外旅游人数下降到了1.78万人,比去年下降了27.22万人,占99.34%。SARS过后海外旅游人数有了较快的回升,在近期有望达到疫前水平。
从历年各月旅游人数的数据可以看出,近年来旅游行业保持着良好势头,中国旅游业的进步是有目共睹的。如果没有SARS,2003年旅游业仍然会保持持续增长。SARS只是一个突发事件,不会对旅游业造成长期的影响,不过在近几个月还有一个恢复的过程。因此认为对2003年今后几月旅游人数的影响因素有:
1、旅游业内在的发展动力。如不断增加旅游业投资,为2008年奥运会做的建设等。
2、SARS流行期间对旅游业的影响。 两种因素的影响特点有:
第一种是旅游业进步的内在本质因素,具有根本性、长期性。第二种是突发事件,是偶然因素,具有暂时性,作用效果看有突出重要性。
分别分析两种力量对2003年旅游业的作用,两种力量的效果是独立影响的,最后把两种效果叠加起来。 (二)模型的建立和求解
1、不考虑SARS对2003年各月旅游人数的预测
从1997—2002年的数据可以看出,各月人数的分布有明显的规律性,见图9 规律一:历年各月都在不断增长,从1999年到2002年4年的平均增幅为27.2万人。增长过程见图10。我们做了每月与历年该月的相关性检验,发现各月与前两年的该月份显著相关。
规律二:每年中各月的大小关系基本不变,表现为图9中的同时上升,同时下降,且各月的增幅比例大致相同。我们也作了12个月中每个月与相邻几个月的相关性检验,发现有一定相关性。如果有足够的数据,我们希望得出这样的相关关系:
F(i,j)=?rkF(i?k,j)???lF(i,j?l)
F(i,j)表示第i年,j月的人数;rk反映历年该月的影响;?l反映相邻月份的影响。但是目前得到的数据太少,不能统计出相关系数。下面从数据的特点预测没有SARS影响的各月旅游人数。
近几年来,旅游业的不断进步,必然吸引越来越多的国外游客。不过各年的增幅有一些波动,其实有点波动是很正常的,因为事物都是在曲折中前进,螺旋式上升的。
从图10中可见,2001到2002年(4点到5点)的增长轨迹与1998-1999-2000年(1-2-3点)的增长轨迹很相似。从1998年开始经过两年快速增长期后,进入一个低速增长期,两个时期的总长约为3年。2001年开始的快速增长,很可能又要进入下一个快速增长期,可推测2003年在2002年的基础上会增长40万人。这40万总增量分配到12个月中,分配比例按表3中的统计比例。得到2003年各月的预计值并加入图9中比较。
11
历年旅游人数分布图403530(万人)25201510500510月15200120022003没SARS预测值 图1999 9到2002各年增幅60(万人)40200-200123456 图10:近几年旅游人数的年增长数 表3:各月份增幅比较 月份 1月 2月 3月 4月 5月 6月 平均增幅 1.025 2.17 1.825 2.25 2.375 2.4 增幅比例 0.041214 0.087254 0.073382 0.09047 0.095497 0.096502 月份 7月 8月 9月 10月 11月 12月 平均增幅 2.05 2.225 2.5 2.025 2.275 1.75 增幅比例 0.082429 0.089465 0.100523 0.081423 0.091476 0.070366 表4:2003年预测值 1月, 2月, 3月, 4月, 5月, 6月 33.3646626.1820332.6997633.0108631.453015.431, , , , , 8 7月, 8月, 9月, 10月, 11月, 12月 35.6219536.0197833.0419825.855329.462, 32.2, , , , 7 2、考虑SARS对旅游人数的影响
用2003年的预测值减去2003年8月以前的实际值,可看出SARS减少的旅游人数,见图11。其中对2003年2月作了处理,2月时北京没有SARS,而且1月、3月的人数比2002年都有增加,符合增长规律。则2月份旅游人数的降低(近一半),由其他突发因素引起。用2002年2月的数据替换。
考虑到SARS的蔓延是每天变化的,把每月的旅游人数离散到每天。把每月人数作平均,作三次样条插值,得到光滑的曲线,曲线上的点反映每天访问人数。
12
差值曲线反映了SARS没天减少的访问人数 。
我们发现每天减少的访问人数有个峰值,这个峰值出现在5月,6月之间,到8月时已经过。说明SARS对5、6月影响最大,到8月已逐渐减弱。关于SARS的指标有多个,如果有个指标能反映对旅游者心理的影响,则SARS对旅游业影响减弱时这个指标也下降。经过对众多指标(死亡率,死亡人数,治愈人数,治愈率,疑似患者数,累计患者数)的筛选,选出能较好地反映对旅游人数的影响的指标——每日现有患者数,并和旅游人数变化曲线画在一起观察,如图12。
图11:曲线比较图
图12
分别选取两条曲线中的两个序列作相关分析,得到每日现有患者数对旅游人数的影响有个时延,这个值为8天,也正反映了两条曲线的峰值相差10天左右。到2003年6月23日,患者还有53个,2个疑似患者,可以确信SARS即将消失,因为已经没有疑似患者,一段时期内多个危险指标没有增长。那么SARS对旅游的影响也即将过去。目前患者数和日访问人数都一直保持恒定的下降速度。可从数据得出日访问人数以0.00976159133019的斜率降低,10月25日以后SARS已没有影响了。
累加SARS对每天访问人数的减少量,得出9月比预计值减少11.5944万,10月比预计值减少2.95万。
13
最后我们综合两种因素后得到预测值见表5,图13。
表5:预测值
1月 15.4 7月 8.8
2月 29.7 8月 16.2 3月 23.5 9月 24.02755 4月 11.6 10月 33.06888 5月 1.78 11月 33.04198 6月 2.61 12月 25.85537 历年旅游人数分布图403530(万人)2520151050020012200246810122003年实际值142003没SARS预测值
图13
以上只分析了旅游业的内在增长因素和SARS的打击作用。在SARS时期没有旅游的人,可能累积到SARS以后旅游,这会使实际旅游人数比预测值大。可见以上只是保守的预测,不过都是很乐观的预测。
五、 短文
传染病模型的研究任重道远
传染病是危及人类身体健康的重要因素之一,长期以来一直受到世界各国的关注。随着社会的发展、卫生设施和医疗水平的改善和人类文明的不断发展,诸如霍乱、天花等曾经肆虐全球的传染病已经得到了有效的控制。但是在世界的某些地区,由于战争对抗、人为失误或病毒变异等因素,还不时出现传染病流行的情况,如82年英国的全国性沙门氏菌病暴发流行,以及最近迅速蔓延全球的SARS等。
目前很多学者对SARS的传播进行了大量的研究。大部分研究采用常微分方程方法进行分析,以个体相互作用的方式研究病毒的传播。最近中国科学院研究
14
生院借鉴流体力学中宏观偏微分方程与微观分子动力学方法相结合的方法,首次建立了传染病传播概率模型。这些模型的出现为疫情的评估和控制提供了科学数据,其重要性可从这次SARS事件中窥见一斑:
1、科学预测,合理分析。
在疫情忽然爆发引起公众普遍关注后,人们最关心的问题有:“SARS疫情将如何发展?何时才能得到有效控制?有没有办法对它进行科学预测?”以上问题均可以通过建立数学模型来解决。作为中国的首都,同时又是重灾区的北京,正是众人关注的焦点。之前已有模型预测北京在实施严格的隔离政策之后,4月底高峰期过后将呈现下降趋势,直至7月将实现零增长,届时累计患者数将达2800多人。事实证明这一预测虽不全中亦不远矣。
2、理智对待,减小恐慌。
在疫情爆发初期,未知病毒的神秘性容易对人造成心理干扰,导致大片人群的恐慌。但根据传染病模型,我们知道当一个传染期内每个病人有效接触易感者的平均人数小于一个人时,病人的比例将逐渐变小,最终实现零增长。也就是说一旦采用严格控制隔离的方法有效地切断了病源和易感人群的联系,将病人的有效接触人数减少到一个以下,必然可以达到有效控制或消除疾病的目的。
3、模拟预测,指导控制措施。
一些模型对政府控制的力度和时机作了相应的模拟,为各地区SARS疫情发展前景预测及政府控制疫情政策的制订给出了科学依据,并提出了一些实用的建议。如:及早采取果断控制措施;即使防范措施晚了一些,出现了高峰,但只要控制严格,降低感染率(即严格进行人群隔离,切断传染途径),前景还是乐观的,因为疫期长度与感染人数峰值并不呈正比例关系等。
4、吸取经验教训,预防疫情反弹。
即使在疫情已接近尾声的今天,对传染病模型的研究仍在继续。在此次SARS战役中,我们获得了宝贵的数据,也吸取了不少的经验教训。通过模型预测,有专家指出:后期绝对不能松懈大意,哪怕是几个病人逃脱控制,也可能造成疫情的反弹,不仅会造成总感染人数的增加,更为严重的是导致疫期的显著延长,由此对我国的国际形象、社会心理、经济增长造成重大损失。
目前对SARS传播的研究主要是以个体相互作用的情况而言,其计算结果的准确性、可靠性将受到限制。随着疫情接近尾声,对SARS的传染病学机理和传播机制越来越了解,可考虑疫情的空间分布和时间发生的概率,结合计算机进行大量模拟,得出更为全面和准确的模型。进一步的工作和更准确结果的得出,将有待于专家学者们进行深入研究。
六、 模型的评价和改进方向
(一)模型的评价:
1、附件1的模型中,未能明确指出自由携带者和已确症患者之间的区别,未能考虑到潜伏期不具备传染力的情况。而在本文中的两个模拟模型,对各个时期的带菌者作了严格的界定,认为只有未受控(未被隔离及入院就诊)的带菌者在病发后的1~2天内才具有传染能力,之后会立即入院治疗。符合实际情况,揭示了传染的主要途径和传播机理。
15
2、改进后的模拟模型对潜伏期和带菌者被发现的时间按照实际规律,作了正态模拟,具有明确的实际意义。
3、模拟模型是完全依照北京市5月10日以前的数据拟合得到的,不需要借鉴其他城市的数据,可以预测首发城市的疫情发展,且结果与实际情况符合的很好。
4、正是由于我们模型逻辑的准确性,如果经过深入调查获取实际数据,再对模型进行一定的修改,就可立即用于建立前期的预测模型,而不需要等待疾病的发展来获得大量数据进行拟和。
(二)理想模型建立的方案及困难
要真正建立能够预测以及能为预防和控制提供可靠,足够的信息的模型。要怎样做,以及困难。
事实上,真正情况要比我们的模型复杂的多,比如,由于初期传播机理不明,医护人员为收到好的保护,很多受到传染。而传播率事实上也是按照某种趋势变化的,他与人们的警惕心理等因素有关。
要建立这样的模型,我们认为:
1、要不断的研究传染病理,完善模型。
2、建立实时反馈机制,不断修改模型参数,以求更准确的预测以后的数据。3、尽量以实际的研究数据建模,力求准确反映事件的实质。这样,当人为修改参数进行模拟预测时,可以为采取什么措施效果最好做出较准确的建议。
4、研究各种措施实际对疾病传播的影响。
要做到这些,也恰恰是建模的难点,因为当涉及到新的传染病的时候,以上数据或参数都必须通过长期的研究或拟合才能得到的。比如措施对疾病传播的影响,是要统计大量实际数据才能获得的。以上也可以看作我们模型的改进方向。
七、 参考文献
[1]范学工 龙云铸,严重急性呼吸综合征——非典型性肺炎讨论,中国现代医学杂志,第13卷第8期,1-4页,2003年4月
[2] 济南市卫生局/机构设置/疾病控制处,SARS传染能力及超级传染
者揭秘,传染性非典型肺炎防治, 2003-6-2
[3] 中国科学院研究生院管理学院“SARS”管理与控制研究小组,对各
种措施的效果分析和结论,研究简报,第1期
16
附件
附件1 确定型模拟程序 function logi(N) Tfault: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; for i=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; for i=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')
17
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
附件2 随机型模拟程序(主程序:carlo)
function P=carlo(M) global a f bb
%脚本程序:clear;M=100;P=carlo(M);sum(P),length(find(P>0.1)) %N是采取措施的天数(4月20日)。M是模拟次数。 %返回的P值是新增患者数的序列。
%N=35,M=500:sum(P)=2831,-----109天 %N=35,M=200:sum(P)=2792,-----111天 %N=30,M=200:sum(P)=1958,-----100天 %N=40,M=200:sum(P)=4487,-----121天
%滞后15天:N=50,M=200:sum(P)=10208,139天
%a是日增患者数。 %f是自由带菌者数。 ?是遗漏率。0
P=zeros(1,200); for j=1:M
f=zeros(1,200); a=zeros(1,200);
18
%f(1)的起点:
f(1:7)=[1,2,3,4,5,6,7]; %初始化自由带菌者数。 bb=0.9;k=3; for i=1:N
a(i+6)=f(i)+a(i+6); func1(k*f(i),i+6); func2(k*f(i),i+7); end
bb=.4;k=1; for i=N:180
a(i+6)=f(i)+a(i+6); func1(k*f(i),i+6); func2(k*f(i),i+7); end
P=P+a; end P=P/M;
nn=length(find(P>0.1)); figure;
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];
plot(37:36+length(AA),AA,'r:');hold; plot(1:nn,P(1:nn),'b--');hold legend('实际','拟和') title('改进后的模拟')
function func2(q,p) global a f bb
%func1---指q人群被感染后当天(感染源就在当天发现)就开始被寻找。(成为隔离的人群)
%b是遗漏率。0
%1-b(1),b(1)*(1-b(2)),b(1)*b(2)*(1-b(3))
bbb=bb^(1/3); ?b为每天的遗漏率,bb为总遗漏率。b是模拟值。 if q<.5 else
n=floor(.5+normrnd(5,2)); if n<1;n=1;end for i=1:3
19
b(i)=normrnd(bbb,.05); if b(i)<0 b(i)=0; elseif b(i)>1 b(i)=1; end end
if n>2
a(p+n)=(1-b(1)*b(2)*b(3))*q; %一部分染病者在被隔离时发病
f(p+n)=b(1)*b(2)*b(3)*q; %另一部分未被隔离,成为自由带菌者
elseif n==2
a(p+n)=(1-b(1)*b(2)*b(3))*q; f(p+n)=b(1)*b(2)*b(3)*q;
func2(b(1)*b(2)*(1-b(3))*q,p+2); %第二天虽被找到,但已造成当天的影响。
elseif n==1
a(p+n)=(1-b(1)*b(2))*q; %已在隔离区中(第零,一天找到的人数和)
a(p+n+1)=b(1)*b(1)*q; %迟一天才找到的(第一天之前(发病之前)未被隔离)。
func1(b(1)*b(2)*q,p+2);%他从发病到入院,有两天的自由感染期。(第二天才被找到)
func2(b(1)*b(2)*q,p+3);
func2(b(1)*(1-b(2))*q,p+2); %第一天虽被找到,但已造成当天的影响。 end end
function func1(q,p) global a f bb
%q是接触患者后受到感染的人数,p是他们被感染的日期。
%func1---指q人群被感染后一天(感染源才被发现)才开始被寻找,(成为隔离的人群).
%b是遗漏率。0
%1-b(1),b(1)*(1-b(2)),b(1)*b(2)*(1-b(3))
bbb=bb^(1/3); ?b为每天的遗漏率,bb为总遗漏率。b是模拟值。 if q<.5 else
n=floor(.5+normrnd(5,2)); if n<1;n=1;end
20
for i=1:3
b(i)=normrnd(bbb,.05); if b(i)<0 b(i)=0; elseif b(i)>1 b(i)=1; end end if n>3
a(p+n)=(1-b(1)*b(2)*b(3))*q; f(p+n)=b(1)*b(2)*b(3)*q; elseif n==3
a(p+n)=(1-b(1)*b(2)*b(3))*q; f(p+n)=b(1)*b(2)*b(3)*q;
func2(b(1)*b(2)*(1-b(3))*q,p+3); %第三天虽被找到,但已造成一天的影响。
elseif n==2
a(p+n)=(1-b(1)*b(2))*q; %已在隔离区中(第一,二天找到的人数和)
a(p+n+1)=b(1)*b(1)*q; %迟一天才找到的(第二天之前(发病之前)未被隔离)。
func1(b(1)*b(2)*q,p+2); %他从发病到入院,有两天的自由感染期。
func2(b(1)*b(2)*q,p+3);
func2(b(1)*(1-b(2))*q,p+2); %第二天虽被找到,但已造成一天的影响。
elseif n==1
a(p+n)=(1-b(1))*q; %分析同上。 a(p+n+1)=b(1)*q; func1(b(1)*q,p+1); func2(b(1)*q,p+2);
func2((1-b(1))*q,p+1); end end
附件3 SARS减少的每日旅游人数 Clear;
Y=zeros(1,160);
for i=1:160;%以-0.00976159133019为斜率 0.00684210526316为切距 Y(i)=-0.00976159133019*i+0.684210526316; end
plot(1:160,Y);%8月15日后SARS对每日旅游人数的减少量,10月25后影响消
21
失。
a=sum(Y(1,1:15))%8月15到8月底 b=sum(Y(1,16:45))%9月受减小量 c=sum(Y(1,46:70))月受减小量。
22