传染病动力学

传染病动力学模型姓名：魏薇薇学号：2009210927院系：数理与信息学院专业：系统理论摘要：本文首先介绍传染病动力学的相关概念,接下来介绍两个基本的传染病动力学模型,最后建立一个传染病动力学的偏微分方程模型,并对模型做一些适当的分析.关键词：传染病动力学;常微分方程;偏微分方程;数学模型Model of Epidemic DynamicsAbstract：This article first introduces the concepts of epidemic dynamics, followed by two basic model of epidemic dynamics, finally it creates a partial differential equations model of epidemic dynamic ,and do some proper analysis to the model. Keywords：Epidemic dynamics;Ordinary differential equations;Partial differential equations;mathematical model前言传染病动力学是对传染病的流行规律进行理论性定量研究的一种重要方法.它是根据种群生长的特性,疾病发生和在种群内传播的规律以及与之有关的社会等因素,建立能反映传染病动力学特性的数学模型,通过对模型动力学性态的定性、定量分析和数值模拟,来显示疾病的发展过程,揭示其流行规律,预测其变化发展趋势,分析疾病流行的原因和关键因素,寻求对其预防和控制的最优策略,为人们防治决策提供理论基础和数量依据.与传统的生物统计学方法相比,动力学方法能更好的从疾病的传播机理方面来反映流行规律,能使人们了解流行过程中的一些全局性态.传染病动力学与生物统计学以及计算机仿真的相互结合、相辅相成,能使人们对疾病流行规律的认识更加深入、全面,能使所建立的理论与防治策略更加可靠和符合实际.1 两个基本的传染病动力学模型在传染病动力学中,长期以来主要使用的数学模型是所谓的“仓室”模型，它的基本思想由Kermack 与McKendrick 创立于1927年,但一直到现在仍然被广泛的使用和不断地发展着.下面我们以他们提出的两个经典的基本模型为例,来阐述建立仓室模型的基本思想和有关基本概念,并显示由模型所能得到的主要结论.1.1 K-M 的SIR 仓室模型所谓SIR 仓室模型就是针对某类传染病将该地区的人群分成以下三类（即三个仓室）：易感者（susceptibles ）类 记为()S t ,表示t 时刻未染病但有可能被该类疾病传染的人数.染病者（infectives ）类 其数量记为()I t ,表示t 时刻已被感染成病人而且具有传染力的人数.移出者（removed ）类 其数量记为()R t ,表示t 时刻已从染病者类移出的人数.设总人口为()N t ,则有()()()()N t S t I t R t =++.K-M 的SIR 模型是一个十分简单粗糙的模型.它的建立基于以下三个基本假设：（1）不考虑人口的出生、死亡、流动等种群动力因素.这意味着考虑一个封闭环境而且假定疾病随时间的变化要比出生、死亡随时间变化显著得多,从而后者可以忽略不计.这样,此环境的总人口始终保持为一个常数,即()N t K ≡,或()()().S t I t R t K ++≡（2）一个病人一旦与易感者接触就必然具有一定的传染力,这里假设t 时刻单位时间内,一个病人能传染的易感者数目与此环境内易感者总数()S t 成正比,比例系数为β,从而在t 时刻单位时间内被所有病人传染的人数（即新病人数）为)()(t I t S β.（3）t 时刻,单位时间内从染病者类移出的人数与病人数量成正比,比例系数为γ,从而单位时间内移出者的数量为()I t γ.显然,γ是单位时间内移出者在病人中所占的比例,称为移出率系数,当不致混淆时也简称为移出率.当移出者中仅包括康复者时,移出率系数又称为恢复率系数或简称为恢复率.在以上三个基本假设下,易感者从患病到移出的过程可用下述框图描述.SIIS I R βγ→→对每一个仓室的人口变化率建立平衡方程式,便得到以下模型：,,.dSSI dt dISI I dt dRI dt ββγγ⎧=-⎪⎪⎪=-⎨⎪⎪=⎪⎩ (1.1.1)下面,我们通过对模型（1.1.1）的分析和解的渐近性态研究来初步显示动力学模型对认识传染病流行规律所起的作用.将（1.1.1）中三个方程两端分别相加,得()0=++dtR I S d ,从而()()()K t R t I t S =++（常数）.由于（1.1.1）中前两个方程中不含R ,故实际上我们只需先讨论前两个方程：()⎪⎪⎩⎪⎪⎨⎧-=-=γββS I dtdI SI dtdS. （1.1.2） 由于0<dtdS,)(t S 单调递减且有下界（为0）,故极限 ()lim t S t S ∞→∞=存在.由（1.1.2）有1dI dSS ρ=-+, γρβ=. (1.1.3)可见,当S ρ=时,I 达到极大值.从而不难在相平面(),S I 上画出系统（1.1.2）的轨线分布图,如图1.1所示.方程（1.1.3）的所有平衡点都在S 轴上,而且0I =为系统（1.1.2）的一条奇线.由图1.1可见,当初始时刻易感者数量()00S S ρ=>时,随着时间增长,染病者数()I t 将先增加达到最大值()I ρ,然后再逐渐减少而最 终消亡.这一现象表明,只要0S ρ>,即011S βγ>,疾病就会流行.Sρ(1.1) 令001S R S βγρ==, (1.1.4)则当01R >时,疾病流行;当01R <时,疾病不会流行,染病者数量()I t 将单调下降而趋向于零.01R =是区分疾病流行与否的阈值.应当指出,（1.1.4）中的1γ表示平均移出时间,也就是平均患病期.事实上,由移出率系数γ的定义可见,若病人数量为n ,则单位时间内移出者的数目为n γ,故经过时间1γ,病人全部移出.要防止疾病流行,必须减少0R 使它小于1,有表达式（1.1.4）可知,这可以通过加强治疗以缩短染病期1γ或采取杀菌等措施以减少疾病的传染力β,或通过隔离措施以减少与患病者可能接触的人数即这里的易感者数0S 来实现.更为有效的方法是通过疫苗接种以使易感者成为免疫者而直接进入移出者类R ,从而减少初始时刻易感者数量0S .设人群中通过接种疫苗成功的比例为()01p p ≤≤,则0S 就变成了()01p S -,从而0R 变小为()0011R p S βγ-=-.要求01R -<,即要求00111p S R γβ>-=-. (1.1.5) 由（1.1.5）式可知,0R 越大,为防止疾病流行所需要接种的人口比例p 就越高.由此可见,对0R 值的估计是十分重要的.由（1.1.4）式可见,要估计0R 的值,难点在于估计β,因为β不仅取决于疾病的种类,而且还依赖于人群所处的社会环境和病人的活动情况.下面介绍一种对0R 的近似估计法.求解方程（1.1.3）,它通过初值()00,S I 的解为()000lnSI I S S S ρ-=--+, （1.1.6） 由于当t →+∞时,()0I t →,()S t S ∞→,代入(1.1.6)式并注意到00S I K +=，得ln0.S K S S ρ∞∞-+= （1.1.7）用数学分析的方法容易验证方程（1.1.7）有且仅有惟一的正实根S ∞.并可解得0,ln ln K S S S γρβ∞∞-==- （1.1.8）0S 与S ∞是可以测定的,例如可以通过血清检查测定.从而可根据（1.1.8）式确定ρ的值,然后由00S R ρ=来确定0R .在测得平均患病期1γ后,也可由（1.1.8）式估算出β.1.2 K-M 的SIS 仓室模型一般来说,通过病毒传播的疾病如流感、麻疹、水痘等康复后对原病毒具有 免疫力,适合用上述SIR 模型描述;通过细菌传播的疾病,如脑炎、淋病等康复后 不具有免疫力,可以被再次感染,1932年Kermack 和Mckendrick 针对这类疾病提 出了康复者不具有免疫力的SIS 模型,疾病的传播机制如下面框图所示：SIIS I S βγ→→这里假设患病者康复后将重新成为易感者,其它假设与SIR 模型相同.此时模型为dSSI I dtdI SI I dt βγβγ⎧=-+⎪⎪⎨⎪=-⎪⎩. （1.2.1） 利用S I K +=,可将方程组(1.2.1)化成方程式()()dS K S S dt βρ=--,.γρβ= （1.2.2）易见,当K ρ≥时,方程（1.2.2）有惟一的平衡点S K =,它是渐近稳定的,即从任一](00,S K ∈出发的解()S t 均单调增加趋向于S K =,从而()I t 将单调减小而趋向于零,说明疾病不会流行.当K ρ<时,方程(1.2.2)有两个正平衡点：,S K S ρ==,S K =不稳定;S ρ=渐近稳定.从任一()00,S k ∈出发的()S t 均随t 的增大而趋向于ρ,从而()1I t ρ→-,这时疾病流行且病人不会消失,最终保持在1ρ-的数量而变成一种地方病.这当然是人们所不希望的.因此,01KR ρ==是区分疾病流行与否,或者是否产生地方病的阈值,当01R <时,疾病逐渐消失;当01R >时,疾病流行而导致地方病产生.2 传染病动力学的偏微分方程模型传染病动力学的常微分方程模型没有考虑到年龄对传染病发展情况的影响.实际上,出生率与自然死亡率如在人口模型中考虑过的那样,应与年龄有关;对传染病本身来说,除极少数疾病（如出血热）外,其发病情况均与年龄有关,有的传染病（如麻疹）在婴儿阶段由于天然免疫力在一段时间内不会发病,同时发病率、治愈率、及死亡率的等也均与年龄有关.因此,必须加入年龄坐标x.此外,病的发展情况通常还和发病时间的长短（病程）有关治愈率和死亡率等均可能和病程有关,因此还需再引入一个病程坐标y.这就使所考虑方程呈现相当复杂的形态.可以考虑下述一些不同的情况：A. 不考虑预防及隔离措施（因而结果偏于“安全”）的情况,或者将预防及隔离措施以某种方式换算为对发病率打一个适当的折扣的情形.a. 病愈后终身免疫的传染病（如麻疹）.b. 病愈后有一段时间免疫力,但不能终身免疫的传染病.c. 病愈后无免疫力,可以立即再感染而重新得病的传染病.B. 考虑预防及隔离措施的情况,或单独考虑其中一种措施的情况.这里又可相应地分为若干情况,不赘述.下面对情况Aa——不计预防及隔离措施,而病愈后为终身免疫的传染病（如麻疹）,建立相应的数学模型.其余情况可类似进行讨论.将全体人口分为三类：Ⅰ.未发病者；Ⅱ.正发病者；Ⅲ.病愈者.三类人之间的相互关系如下图：以下记t为时间,x为年龄,而y为病程.设A 为人的最大寿命,B 为最大病程（B ≤A ）.于是求解区域应为{}B y A x t ≤≤≤≤≥0,0,0.以 ()x t p ,1记第Ⅰ类人在 t 时刻按年龄x 的分布密度函数,()x t p ,3 记第Ⅲ 类人在t 时刻按年龄x 的分布密度函数,()y x t p ,,2记第Ⅱ类人在t 时刻按年龄x 及病程y 的分布密度 函数,于是,在时刻t ,年龄在[]dx x x +,中第Ⅰ类人数为()x t p ,1dx ,年龄在[]dx x x +,中第Ⅲ类人数为()x t p ,3 dx ,年龄在[]dx x x +,、病程在[]dy y y +,中第Ⅱ类人数为()y x t p ,,2dx dy ,其中()x t p ,1及()x t p ,3的定义域为{}A x t ≤≤≥0,0,而()y x t p ,,2的定义域为{}B y A x t ≤≤≤≤≥0,0,0.由于年龄为x的人其病程x y ≤,故当x y ≥时恒有()0,,2≡y x t p决定了函数()x t p ,1,()x t p ,3,及()y x t p ,,2,就决定了此传染病的动力学特征.下面推导它们应满足的方程.注意到对任何人来说,时间增量=年龄增量=病程增量, 我们有在dt t +时刻,年龄在[]dx x x +,中的第Ⅰ类人数()dx x dt t p ,1+应等于在t 时刻年龄在[]dt dx x dt x -+-,中的第Ⅰ类人数()dx dt x t p -,1减去在[]dt t t +,中年龄在[]dt dx x dt x -+-,中的自然死亡数()()dxdt dt x t p dt x d --,1及发病数∧α()dt x t p -,1dx dt .由此可得到()x t p ,1满足()()()()111,,,p p t x t x d x p t x t x α∧∂∂⎛⎫+=-+ ⎪∂∂⎝⎭,（2.1） 这儿()x d 为自然死亡率，∧α为发病率.考虑到传染病的特点，在[]dt t t +,中年龄在[]dx x x +,中的发病人数与人数()x t p ,1dx 及时间dt 均应成正比,同时还和第Ⅱ类人的总数()()⎰⎰=A Bdxdy y x t p t p 0022,,（2.2） 成正比，故()()()()⎰⎰==∧A Bdxdy y x t P x t p x 0022,,ααα,从而（2.1）式可写为()()()()()x t p t p x x d x p tp ,1211α+-=∂∂+∂∂, （2.3） 而()t p 2由（2.2）式定义.同理,对()y x t p ,,2我们有在dt t +时刻,年龄在[]dx x x +,、病程在[]dy y y +,中的第Ⅱ类人数()dxdy y x dt t p ,,2+应等于在t 时刻,年龄在[]dt dx x dt x -+-,、病程在[]dt dy y dt y -+-,中的第Ⅱ类人数()dxdy dt y dt x t p --,,2减去在[]dt t t +,中年龄在[]dt dx x dt x -+-,、病程在[]dt dy y dt y -+-,中的第Ⅱ类人的自然死亡数)(dt x d -()dt y dt x t p --,,2dx dy dt 、传染病死亡数()dt y dt x d ---,()dt y dt x t p --,,2dx dydt 及治愈数()dt y dt x --,β()dt y dt x t p --,,2dx dy dt .注意到()()dt y dt x t p y x dt t p ---+,,,,22=()()()y x t p y x dt t p ,,,,22-++()()()y dt x t p y x t p ,,,,22--+()()()dt y dt x t p y dt x t p ----,,,,22=()()()222,,,,,,,p p pt x y t x y t x y dt t x y ⎛⎫∂∂∂++ ⎪∂∂∂⎝⎭故得()y x t p ,,2应满足的方程为()()()y x t yp y x t x p y x t t p ,,,,,,222∂∂+∂∂+∂∂=()()()()y x t p y x y x d x d ,,,,2⎪⎭⎫ ⎝⎛++--β （2.4）同理，我们有在dt t +时刻、年龄在[]dx x x +,中的第Ⅲ类人数()dx x dt t p ,3+等于在t 时刻、年龄在[]dt dx x dt x -+-,中的第Ⅲ类人数()dx dt x t p -,3减去在[]dt t t +,中年龄在[]dt dx x dt x -+-,中的第Ⅲ类人的自然死亡数()()d x d t dt x t p dt x d --,3加上在[]dt t t +,中年龄在[]dt dx x dt x -+-,中的第Ⅱ类人的治愈数为()()dxdt dy y x t p y dt x B ⎪⎪⎭⎫⎝⎛-⎰02,,,β. 由此我们得到()()()33320,,,Bp p d x p x y p t x y dy t x β∂∂+=-+∂∂⎰（2.5） 下面看初始条件及边界条件. 初始条件为:0=t ()x p p 011=,()y x p p ,022=,()x p p 033=. （2.6）边界条件:由于新生婴儿进入第Ⅰ类状态，且从o x =开始，故出生的婴儿数将给出在0=x 的边界条件.设出生率为()x b ，并设最低生育年龄为()A a <,我们得到在时段[]dt t t +,中出生的婴儿总数()()()()dxdt x t p dy y x t p x t p x b Aa B⎰⎰⎪⎪⎭⎫ ⎝⎛++0321,,,,应等于在时刻dt t +、年龄区间在[]dt ,0中第Ⅰ类人数()()()dt t p dt dt t p 0,0,11=+, 故有边界条件:0=x ()0,1t p =()()()()ξξηηξξξd t p d t p t p b BAa ⎪⎪⎭⎫⎝⎛++⎰⎰0321,,,,, ()0.t ≥（2.7） 此外，第Ⅱ类及第Ⅲ类人中不包括新生婴儿,故应有边界条件:0=x ()0,0,2=y t p ()B y t ≤≤≥0,0, （2.8） :0=x ()00,3=t p ()0≥t . （2.9） 又由于第Ⅰ类人中的发病者进入第Ⅱ类人，病程从0=y 开始,故第Ⅰ类人的发病数应给出0=y 处的边界条件.我们有在[]dt t t +,中年龄在[]dx x x +,中的第Ⅰ类人的发病数()dxdt x t p ,1∧α应等于在t 时刻,年龄在[]dt dx x dt x -+-,、病程在[]dt ,0中的第Ⅱ类人数()dxdt dt x t p 0,,2-,故得如下的边界条件 :0=y ()()x t p x t p ,0,,12∧=α ()A x t ≤≤≥0,0. （2.10） 其中()()t p x 2αα=∧,而()t p 2由（2.2）式定义.这样就得到定解问题（2.3）-(2.10),其中()t p 2由（2.2）式定义.在对已知的资料加以适当的假设（包括相容性条件）后,我们要求该问题的解()x t p p ,11=,()y x t p p ,,22=,()x t p p ,33=,使在区域{}B y A x t ≤≤≤≤≥0,0,0上连续.可以看到这个问题有如下一些特点：(ⅰ)有三个未知函数，其中()x t p ,1及()x t p ,3具有两个自变数,而另一个未知函数()y x t p ,,2则具有三个自变数.它们的方程及边界条件互相耦合在一起.(ⅱ)(2.3)-（2.5）均为主部为常系数的一阶偏微分方程,但除(2.4)是关于()y x t p ,,2本身（无耦合）的普通的一阶线性偏微分方程外,关于()x t p ,3的方程(2.5)中包含2p 对y 的积分,因而是线性、非局部的方程,而关于()x t p ,1的方程（2.3）由于包含()t p 2,不仅具非局部的形式,而且是非线性的.（ⅲ）在0=x 处对1p 的边界条件具有非局部的积分泛函的形式,但还是线性的;而在0=y 处的边界条件不仅是非局部形式,而且是非线性的.总之,这是一个一阶双曲型方程组的非局部、非线性混合初边值问题,而且未知函数具有不同个数的自变数.对这类方程组的定解问题尚有大量问题（如解的整体存在性、解的性质等）有待进一步研究讨论.结束语用数学方法来考察传染病的理论,对它的发病机理、动态过程和发展趋势进行研究,已逐渐成为一个活的研究领域.本文首先介绍了两个经典的传染病动力学模型,然后引入多个变量从偏微分方程的角度来考察传染病的流行规律.从而使所建立的模型与实际更加符合,也能更好的研究传染病的流行规律.。