结核病传播的数学模型(代胜吉焦胜博朱益飞).doc

结核病传播的数学模型代胜吉 焦胜博 朱益飞摘要本文建立结核病传染的改进K-M仓室模型，并对广州市结核病传播情况进行了分析和预测；讨论了不同预防和控制措施对结核病传染的影响；最后考虑了人口流动对模型的影响，将模型进一步完善。首先，将人群分为7类并确定了他们之间的关系，建立了改进型K-M仓室模型。列出了传染病动力学微分方程组和差分方程，计算并绘出了广州市地区结核病传染相关指标 （至每一季度季度累计染病人数）、（至每一季度累计无法治愈人数）的理论曲线，将其与实际曲形进行了比较，得到了比较满意的结果。根据2008年广州市第一季度初始值预测了2008至2012年广州市结核病累计染病人数（），得出结论：到2012年，广州累计结核染病人数将达168940例，累计不可治愈人数将达173660例。接着，分析了政府采取的措施对结核病传染的影响。根据各措施对模型造成的不同影响将措施分成四类并分别分析其对结核病传播的影响，提出了改进意见：（1）着重建立和完善结核病监控体系，切实实行好各项措施；（2）扩大对结核病治疗的范围，增加用于治疗结核病的经费，提高对耐药性结核病的治疗力度；（3）艾滋病的检测和治疗对结核病流行趋势的抑制效果不能立刻体现，因此在采取各项预防和控制结核病措施的同时还应常备不懈地坚持对艾滋病的检测和治疗；（4）在控制艾滋病增长的同时，扩大艾滋病免费治疗的范围，对艾滋病和结核病交叉感染的治疗实施政策优惠。最后，将流动人口这一因素考虑到模型之中，对问题一中建立的模型进行进一步改进。经计算得出结论：（1）流动人口对结核病传染的影响是非常大的，人口的流动会极大促进结核病的传播。在考虑外来流动人口的情况下，结核病流行程度比无外来流动人口的情况严重，2005年末前一种情况的总染病人数（145660人）已是后一种情况（43362人）的约3.36倍，且这种倍数关系随着时间推移有增大的趋势。（2）人口流动还会影响正常的结核病预防控制措施的实施效果。关键词 K-M仓室改进模型 结核病 预防控制 流动人口一、问题的重述结核病是一种古老的经呼吸道传播的慢性传染病，在历史上流行猖獗，被称为白色瘟疫。仅从1882年著名的德国科学家郭霍氏发现结核菌以来，至少有约2亿人被结核病夺去了生命。尽管上世纪50年代有效的抗结核药物问世以后，形势发生了根本的改变，但是进入八十年代以后，由于各国政府的忽视、移民、多耐药结核病和艾滋病等因素的影响，在全球，特别是在发达国家又重新出现了上升的趋势。我国是全球22个结核病高负担国家之一，结核病人数位居世界第二位，仅次于印度。我国卫生部于1979年至2000年在全国先后开展了四次结核病调查，结核病发病数居法定传染病之首。2010年我国进行第五次结核病流行病调查。2000年结核病流行病调查结果分析：我国结核病疫情还相当严重，部分地区有蔓延趋势， 表现为结核菌感染人数多、现患肺结核病人多、结核病死亡人数多、农村结核病人多和传染性肺结核病疫情居高不下。结核病有着独特的传播规律。它的流行与其它的传染病非常不同。它的病情发展很慢、在疾病的潜伏期和活化期采取完善的措施完全可以根治，但如果治疗不当或治疗不彻底而发展为耐药性结核就会成为不治之症，而且结核病病愈后不会获得免疫力。目前，在全球范围内，由于耐药结核菌的产生与扩展，结核菌与人体免疫缺陷病毒(HIV)的双重感染以及许多国家结核病控制规则的不完善，使得全球结核病商情呈明显上升趋势。在全球所有传染性疾病中，结核病已成为成年人的首要死因。结核病流行的另一在特点就是其潜伏期很长，而且因人而异，有的感染者也可能终生都不会发病。查询相关资料解决以下四个问题： 1、建立你们认为能反应我国或某一地区(不限中国国内)的结核病传播规律的传染病动力学模型，并对其结核病流行趋势进行中长期预测；2、通过模型分析，对现有的预防和控制措施进行评价，并分析不同的预防和控制措施对结核病传播的作用效果，也可提出自己的最佳预防控制策略；3、尝试将人口流动特点(如：我国的农民工在城乡间的流动，或不同地区国家之间的人口迁移)等因素对结核病传播和发作的影响予以定量化描述，并引入到所建模型之中，寻求对结核病传播和预防控制的影响。4、 给当地报刊写一篇通俗短文，说明对结核病进行预防控制的重要性。二、模型假设1、非耐药性结核病治疗周期极短，且治疗后无传染性；2、认为未接种疫苗人群即为易感染人群；3、若艾滋病与结核病交叉感染则无法治愈；4、不考虑人口流动对艾滋病增长的影响；5、不考虑所研究地区的出生死亡率；6、每年所研究地区的净流入人数为定值；7、选定地区的人口年龄结构在所研究时间段里基本不变；8、只有儿童必须依政策进行结核病免疫接种。三、符号说明： 易感人群数（未接种疫苗）： HIV/AIDS人数： 未及时治疗或治疗不当的发病人数： 发病人群中治疗治愈人数： 不可治愈人数（耐药性、HIV交叉感染）： 所研究地区总人数： 未及时治疗及治疗不当者所占总感染细菌者的比率： 每个季度每个染病者平均作次能引起传染的接触： 人口增加影响系数： 易感人群受传染的传染率： HIV感染者中感染肺结核细菌的比率： 每季度HIV的死亡人数： 每季度肺结核的死亡人数： 单位季度HIV增长率： 每季度净流入广州的外来人数： 每季度净流入广州的外来人数的患病比率四、问题1 初步模型建立和分析4.1初步模型的建立在建立初步模型时，考虑到结核病具有不同于一般传染病的特点：（1）未及时治疗或治疗不当的结核病患者会发展为耐药性结核无法治愈；（2）HIV患者易感染结核病且无法治愈；基于以上特点，我们选择对传统的K-M“仓室”( compartment) 模型进行分析并及加以改进。4.1.1人群的分类:（1）不易感人群（接种疫苗不发病人群）（2）易感人群 （3）HIV/AIDS人群（4）未治疗和治疗不当者（5）无法治愈者（耐药性、AIDS交叉感染等）（6）治愈者（7）死亡者以上几类人群的关系如改进后的仓室模型框图1所示：HIV/AIDS人群未治疗和治疗不当无法治愈（耐药性、AIDS交叉感染等）易感人群死亡治愈图1 改进后仓室模型框图4.1.2 初步模型的建立根据图1所示模型框图建立结核病动力学微分方程组如下：（4）（1）（2）（3）及艾滋病人数随季度增长的差分方程：（5） 4.2模型的求解（1）确定研究对象考虑到作为改革开放后首批发展迅速的沿海城市，广州经济发展十分迅速，但在结核病传染防治等方面也表现出了其非常脆弱的缺点：人口复杂、成分多且结核病病毒携带者病例年年增加。为此，选取具有代表性的广州市为研究对象。（2）求解微分、差分方程求解所建立的微分方程组其解析解是十分困难，我们根据搜集的数据，运用MATLAB软件给出了微分方程组的数值解，并推导了差分方程的解析解。4.2.1模型初值的确定查01年第一季度广州未接种疫苗人群数据知；查01年第一季度广州登记肺结核病治愈人数；4.2.2推导和计算（1）求出差分方程（5）的解析解并代入微分方程组中。根据参考资料【】查得；（2）利用MATLAB工具箱求解微分方程组（）得到、的数值解；4.3模型结果分析（2001-2006年）4.3.1对结核病流行趋势评断指标的确定不失一般性，设定以下两指标对广州市结核病流行趋势进行分析和预测：（1）：至该季度总染病人数；（2）：至该季度无法治愈人数；4.3.2对2001-2006年广州市结核病流行趋势分析的计算结果经计算，得到结核病流行趋势按季度变化结果如表1所示：表1 模型数值计算解年季度（人）（人）200119834 10024 211816 12036 313706 13960 414716 15013 2002115023 15370 215815 16223 317446 17927 419065 19634 2003120223 20896 221377 22175 322965 23915 424869 26000 2004126738 28087 228514 30124 330483 32406 432778 35077 2005135247 37998 237777 41069 340442 44383 443362 48083 2006146559 52215 249960 56739 353556 61682 457426 67168 为了便于观察，我们将模型计算结果和根据国家公共信息网公布数据得到的流行趋势绘制成图2进行比对如下：图2 计算数据与实际染病数对比图对比模型计算结果和实际结果可见，在第1至第20季度计算结果和实际结果基本吻合，第20季度以后误差开始显著增大，即该模型可以较好的预测五年之内的结核病流行趋势。4.4未来5年结核病传染预测4.4.1预测时间长度的确定从1953年第一个国家五年计划至今，我国已经完成了11个五年计划。国家对一个计划的周期进行了科学的选择，计划包含了医药卫生在内的许许多多重要领域，对我国的快速发展起到了巨大的促进作用。我们根据国家五年计划这一规划，将广州市的结核病传染中长期预测定位为5年，五年的医疗卫生计划既符合国家的发展计划，针对结核病传染也可以科学地采取相应的措施，对于政府的决策也用很大的指导意义。4.4.2未来五年结核病传染预测利用广州2008年第一季度数据作为初值，带入初步模型中，计算结果列于表2：表2 模型预测值（2008-2012）年季度2008145318455082481314833635133651464455087552122009159143595352634196400136800468498472938733812010178229788512838998476438998290987496506976042011110350010479021110101126103119060121010412769012999020121136940139660214687015014031575201614604168940173660流行趋势预测曲线如图3所示图3 2008-2012总染病人数预测根据初步模型计算结果和预测曲线可见，如果不采取有效措施，随着时间推移，结核病总染病人数将会迅速增长。4 .4模型的评价在考虑艾滋病的影响和结核病存在耐药性两因素的基础上，给出了广州市结核病发病及传染的模型，并做出了从2008年到2012年的流行趋势预测，给出了结核病总染病人数和无法治愈人数随时间推移的发展趋势和具体计算结果，对政府及医疗机构的科学决策提供了依据。该模型没有考虑人口出生死亡对人口总数的影响、医疗水平提高对治愈人数的影响等因素，且认为艾滋病按常指数增长，这些都影响了使用该模型中长期预测的准确程度。五、问题2 预防和控制措施对结核病传播的影响结核病的传播对人们的身体健康、社会安定、经济的发展等方面有着深远的影响，采取科学的预防和控制措施对结核病的传播具有巨大的遏制作用。采取不同的预防和控制措施对结核病的传播具有不同程度的遏制，高效的措施将发挥其不可替代的独特作用，对控制结核病传染起到“事半功倍”的效果。我国目前采取DOTS策略来预防和控制结核病的传播，主要包括以下措施：（1）政府对结核病规划的承诺；（2）通过痰涂片镜检发现病人；（3）在正确的管理下，给予标准的的短程化疗；（4）建立正规的药物供应系统；（5）建立对规划执行的监督，评价系统；问题1中建立了分析预测结核病传播规律的初步模型，下面以初步模型为依据，对现有的预防和控制措施进行评价，并给出改进意见：5.1措施评价指标的确立为便于评价所采取的预防控制措施对结核病流行情况的影响，仍将季度总染病人数和季度无法治愈人数作为评价指标，并与未采取措施的情况进行对比。5.2现有预防和控制措施对结核病传播的影响预防和控制措施对于初步模型的影响，主要体现在对初步模型中的参数、未及时治疗及治疗不当者所占总感染细菌者的比率、单位季度HIV增长率、HIV感染者中感染肺结核细菌的比率值的变化。故分别这四种情况加以分析：5.2.1 改变的措施及其效果经过调查，现阶段改变的措施主要有以下4项：（1）通过痰涂片镜检发现和治疗染病人群以切断传染源以降低有效接触次数；（2）加强对易感人群（艾滋、弱体制、特定年龄段）的保护、定期体检以降低易感人群受传染的传染率；（3）加强宣传力度，督促养成良好卫生习惯以降低和；（4）建立特定人群（学生、工人、军人等）定时检疫和接种疫苗的制度以降低易感人群数量初值；为研究减小对结核病流行情况的影响效果，我们将该参数值缩小后代入初步模型中计算2001至2006年总染病人数，并将其与参数值缩小前的总染病人数对应的绘制在图4中：图4 采取减小措施前后总染病人数对比图从图4中可见，与未采取措施的情况相比较，采取减少的策略，可以将总染病人数控制在约左右，明显的抑制总染病人数急剧增加的趋势。5.2.2改变的措施及其效果改变的措施主要如下：（1）对结核病人采取免费治疗的政策以降低（2）对医院进行正规化管理和监督，提高治疗效果以降低在验证降低对结核病流行趋势影响时，我们将该值缩小后代入初步模型中计算2001至2006年总染病人数，并将其与参数值缩小前的总染病人数对应的绘制在图5中：图5 采取降低措施后总染病人数对比图从图5中可以发现，5年（20个季度）后总染病人数由左右降至左右，可见值即使有微小的降低，也会较明显的抑制结核病流行的趋势。5.2.3改变的措施及其效果改变的措施主要有：（1）免费进行艾滋病检测（2）对艾滋病进行免费治疗（3）加强对艾滋病危害的宣传以上三项措施都会降低值，故我们将减小15%后代入初步模型中计算2001至2006年总染病人数，并将其与参数值缩小前的总染病人数对应的绘制在图6中：图6 采取降低措施后总染病人数对比图研究图6知，减小可以在一定程度上减缓结核病流行趋势，但效果需要约15个季度后才较为明显，故对艾滋病的控制应常抓不懈。5.2.4改变的措施及其效果可以改变的措施主要如下：（1）加强对艾滋病人的教育和指导（2）加强对艾滋病人的健康保护以上两项措施都会降低值，故我们将减小70%后代入初步模型中计算2001至2006年总染病人数，并将其与参数值缩小前的总染病人数对应的绘制在图7中：图7 采取降低措施后总染病人数对比图由图7可见，值即使减小较大，对结核病流行的趋势影响仍然较小。由式（）知，这主要是由于艾滋病人数增长较快、基数较大抵消了值减小的效果造成的。5.4科学预防和控制措施的几点建议根据以上分析可以发现，、值的减小对抑制结核病流行趋势有较好的效果，为此，我们在现有预防和控制结核病措施的基础上提出了可以改进方案：（1）考虑到减小对结核病流行趋势的抑制效果最为明显，故应着重建立和完善结核病监控体系的建立，切实实行好5.2.1所列4条措施；（2）的减小程度主要依赖于政府政策对于传染性结核病人的扶住程度。考虑到现阶段政府对于结核病人的免费治疗还仅限于首次得病治疗和非耐药性结核病治疗，故建议政府加强政策支持，扩大对结核病治疗的范围，增加用于治疗结核病的经费，提高对耐药性结核病的治疗力度。这将对降低值起到较好的效果；（3）值对结核病流行趋势的抑制效果不能立刻体现，因此在采取各项预防和控制结核病措施的同时还应长期坚持对艾滋病的检测和治疗；（4）在控制艾滋病增长的同时，扩大艾滋病免费治疗的范围，对艾滋病和结核病交叉感染的治疗实施政策优惠。六、问题3 模型的改进6.1初步建立的模型的不足初步模型将研究对象作为一个封闭的系统进行考虑，没有考虑到外来流动人口涌入的情况，而实际上流动人口尤其是农民工的涌入会对模型产生诸多影响：（1）流动人口无法接受规范治疗，从而增大染病总人数值根据有关资料，农民工已成为新发肺结核病的主要人群，因为农民工群体缺乏结核病防治知识，自我防范意识较弱，他们患病后难以纳入规范化的治疗过程，治愈率较低，一旦脱离治疗就会形成新的传染源在农民工内部和广州本地市民之间传播。这部分流动人口涌入城市，将会对染病人群的值产生较大影响；（2）流动人口流动性强，会增大有效接触次数的值农民工群体一方面由于不是本地常住居民，具有极大的流动性；另一方面存在工作环境的通风、卫生条件普遍不高的缺点，这些必然带来有效接触次数的增加。 而在初步模型中，中没有考虑外来人口，且认为是一个常数，如果考虑外来流动人口，该模型就存在不足了。 6.2模型改进的几点建议 根据以上分析，我们对初步模型进行如下改进：（1）在初步模型中的项中引入农民工中的染病人数，其中为每季度净流入广州的外来人数的患病比率，为每季度净流入广州的外来人数；（2）根据参考资料，引入将初步模型中的常数改为与总人口相关变量，其中为人口增加影响系数。 6.3考虑人口流动的结核病传染模型 6.3.1人口流动结核病传染模型建立改进后的初步模型如下图8所示：图8 改进后的初步模型框图根据改进后的初步模型框图8，列写改进后的初步模型结核病动力学微分方程组如下：（9）（8）（7）（6）以及艾滋病人数随季度增长的差分方程（10）6.3.2模型求解初值的确定（1）查01年相关数据，在第一问的基础上考虑外来人口染病数量，确定01年第一季度；其中假设为0.5。（2）查01年第一季度登记肺结核病治愈人数；推导和计算（1）依初步模型的方法求出差分方程（10）的解析解并带入微分方程组；（2）利用matlab工具箱求解微分方程组得到、的数值解；6.3.3模型结果分析对结核病传播趋势评断指标的确定依然采用初步模型中设定的两个指标、来分析广州市人口流动特点对结核病传播和发作的影响。对2001-2006年广州市结核病流行趋势分析的对比计算结果我们将计算得到的考虑人口流动时结核病流行趋势按季度变化结果和经初步模型计算的不考虑人口流动时的计算结果同时列于表3，并将对比计算结果在坐标图9中表示出来：表3 改进后初步模型的计算结果及比较年季度I1+I2I3考虑流动染病人数不考虑流动染病人数考虑流动染病人数不考虑流动染病人数2001112584 9834 10949 10024 217671 11816 16937 12036 322463 13706 22242 13960 427164 14716 27205 15013 2002131308 15023 31012 15370 235755 15815 35136 16223 341062 17446 40523 17927 446858 19065 46521 19634 2003152681 20223 52320 20896 258609 21377 58035 22175 365053 22965 64362 23915 472153 24869 71522 26000 2004179649 26738 79054 28087 287420 28514 86725 30124 395642 30483 94818 32406 4104490 32778 103610 35077 20051113940 35247 113060 37998 2123930 37777 123000 41069 3134460 40442 133420 44383 4145660 43362 144520 48083 20061157610 46559 156410 52215 2170310 49960 169040 56739 3183750 53556 182380 61682 4198020 57426 196550 67168 图9 模型改进前后总染病数对比图 由计算结果分析人口流动特点对结核病传播的影响由表3和图9可见，在考虑外来流动人口的情况下，结核病流行程度比无外来流动人口的情况严重，如在第20个季度（第5年）时，前一种情况的总染病人数（145660人）已是后一种情况（43362人）的约3.36倍，且这种倍数关系随着时间推移有增大的趋势。根据模型2和计算数据说明，人口的流动会促进结核病的传播。6.4人口流动特点对结核病预防控制的影响在问题二中，我们分析并证明了现有的结核病预防控制措施的效果主要通过减小、来抑制疾病的流行，而根据式1和式6可知人口流动特点的引入，会增加和的值，从而削弱正常的预防控制措施。 6.5改进后模型的评价改进后的模型用有效接触率和总染病人数的变化来表示人口流动特点对模型的影响，计算了人口流动对结核病传播的影响，并分析了其对现有预防控制措施的影响，较问题一中建立的初步模型有更好的适应性，为进一步制定预防和控制措施提供较好的参考。给报刊的短文科学预防控制结核病传染的重要性结核病传染模型启示 自人类出现在地球上，就一直受病毒的侵扰，传染病肆虐成狂。结核病就是一个非常典型的例子，它是一种古老的经呼吸道传播的慢性传染病，在历史上流行猖獗，被称为白色瘟疫。尽管上世纪50年代有效的抗结核药物问世以后，形势发生了根本的改变，但是进入八十年代以后，由于各国政府的忽视、移民、多耐药结核病和艾滋病等因素的影响，在全球，特别是在发达国家又重新出现了上升的趋势，为此，2010年我国再一次进行了结核病普查。针对我国2000年结核病流行病调查结果分析：我国结核病疫情还相当严重，部分地区有蔓延趋势， 表现为结核菌感染人数多、现患肺结核病人多、结核病死亡人数多、农村结核病人多和传染性肺结核病疫情居高不下等。它对我国的经济、社会安定等方面起到了很大的冲击作用，严重影响了我国的国民身体健康，结核病已成为了我国人口非自然死亡的首要因素！传染病对于人们来说一直是一个难题。对于一种传染病，医学工作者会不断深入的探索与研究，然后研制出相关的疫苗或是制定相关的措施，最大程度地减小疾病的危害。当然，在研究传染病本身的同时，我们还迫切需要了解传染病的传染形式、与其相关的措施对病况的影响及未来的传染趋势进行预测等等，以防它对国家人民的健康、社会经济的发展产生更加深远而不可逆的影响。而这些就要依靠数学工具进行处理解决：数学建模在基于现实情况的基础之上对问题就做出了详细的解答，指导人们实践，对传染病传染做出了科学的解释与预测，有效遏制了传染病的进一步发展。我国结核病发展现状令人担忧，防治任务艰巨，结核病的预防控制压力空前之大。从相关的数据上可以很明显的看出来近三十年来我国的结核病患者是逐年增高，由于部分地方政府承诺落实不足、结核病防治经费投入仍缺口较大；流动人口的不断增多、青壮年体质的下降、结核菌/艾滋病毒双重感染和耐药结核病的挑战等因素的影响，这种趋势也是相当明显的，所以，科学的分析与预测与防控就显得更为重要！假若没有了科学的理解，对重要的措施如外来人口的医学检测、当地医疗水平没有进行加强等，没有对传染趋势做出科学的预测，那这些因素将直接影响我国的结核病传染形势，我们可以想象一下现在结核病会对我们的国家社会造成多么大的影响！结核病的预防控制迫在眉睫！对于结核病，有几个主要因素造成了传染趋势上升？又是什么因素造成的结果最严重？政府和医疗机构应各采取多大的力度来处理结核病？这些问题的结果对于结核病的防治是相当重要的也是非常必要的。如果能够科学的回答这些问题，那么将对我国今后的防治结核病工作产生巨大的推进作用！由此可见像我国这样的结核病传染大国，科学的预防控制非常重要。只有我们在平时就做好防范，对于没有根除的疫情，当其重新抬头时，我们就可以从容面对，战胜疫情，并最终将现在还没有根除的传染病根除，为我们社会、国家、世界做出应有的贡献！结核病给我们的正常生产生活造成了很大的影响，2010年全国结核病普查再一次给我们敲响了警钟。我们应不断加强结核病的预防与控制，并以此来指导生产实践，因为它关系着我们每个人的健康！七、参考文献1 王正林、刘明，精通MATLAB7，北京：电子工业出版社，2008。2 余洋、胡志兴，一类具有常数移民且带隔离项传染病模型的分析，石家庄学院学报，第8卷第6期：43-48，2006。3 姜启源、谢金星、叶俊，数学模型（第三版），北京：高等教育出版社，2003。4 杨秀香，接触率依赖于总人数的传染病动力学模型，渭南师范学院学报，第22卷第2期：17-19，2007。5 公共卫生数中心，20012008年广州市结核病监测数据，/CdcDataShare/inceptAction.do?method=doShow，2010年。6 新浪网，我国结核病有六多，，2010。7 赵介民，肺结核的基本知识，北京：科学技术出版社，2009附录1.1 解初步模型中的微分方程function dx=Q1equation(t,x)dx=1;1;dx(1)=1.01*x(1)-x(1)+1.01*x(2);dx(2)=0.06*x(2)-x(2)+0.06*x(1);endt=1:0.1:24;n=length(t);x0=9834 2330;opt1=odeset(RelTol,1e100,AbsTol,1e100);%解微分方程组在此t,x=ode23(Q1equation,t,x