基于Matlab的传染病动力学模型仿真平台

基于Matlab 的传染病动力学模型仿真平台
Simulation Platform of Epidemic Dynamics Model Based on MATLAB 摘要：开发了基于Matlab 的传染病动力学模型仿真平台，通过对传染病动力学模型进行动态仿真，可以对传染病动力学模型的变化进行观察和分析，同时在该仿真平台上，采用时滞微分方程、脉冲微分方程等数值算法实现对传染病模型进行数值模拟，是一个十分实用、方便的仿真操作平台。

关键字：传染病动力学模型；数值仿真；Matlab ；时滞微分方程
Abstract ：A simulation platform of epidemic dynamics model is developed by using Matlab.The simulation platform can be used in dynamics simulation of epidemic dynamics model,and the simulation could be used in results analysis.Based on the platform,Delay differential equations and impulsive differential equations of numerical algorithm can be used to numerical simulation of epidemic model and the multivariable control system simulation.
Keywords ：Epidemic Dynamics Model ，Numerical Simulation ，Matlab ，Delay Differential Equations
1引言
近年来，作为传染病研究的手段之一，利用计算机对传染病动力学模型进行数值仿真越来越受到人们的重视。

诸如MATLAB 中ODE45、DDE23等程序包，被人们普遍使用于传染病动力学模型的仿真中。

近年来随着研究工作的深入，大量新的模型也逐渐受到人们的重视，如：时滞微分传染病模型；脉冲传染病模型；常微分、偏微分混合的传染病模型等。

由于ODE45、DDE23等程序包不是针对传染病动力学模型所开发，无法解决以上这些模型的仿真问题，这些都给相关研究工作造成了一定的困难。

本文利用MATLAB 提供的图形化用户界面（GUI ），结合时滞微分方程、脉冲微分方程等数值算法，并考虑传染病动力学模型的实际研究情况，开发了一套简单、实用的传染病动力学模型数值仿真平台。

2传染病动力学模型的建立
从模型的数学结构来看，传染病动力学模型分为常微分模型、时滞微分模型、脉冲微分模型和偏微分模型等多种形式。

以下以脉冲接种作用下的时滞传染病动力学模型为例，介绍模型的建立方法。

“时滞”可以反映传染病的潜伏期，患者对疾病的感染期和恢复者对疾病的免疫期等实际现象，因此使用“时滞”模型更贴近实际。

如Cooke 等人将时滞因素引入到SEIRS 传染病模型中，用时滞项来反映传染病的潜伏期，建立了如图1
所示的仓室框图。

图1SEIRS 模型的仓室框图
在此模型中，将传染地区的人群分为四类：用S(t)，E(t)，I(t)，R(t)分别表示t 时刻易感者、在潜伏期的感染者、染病者和移出者的数量。

箭头所指方向可以清楚的显示出各类人群流动的情况，τ>0是模型的时滞项，代表疾病在人群中的潜伏期，r >0表示感染者被治愈后返回到易感人群中的速率，β>0是传染率系数，δ为感染者被治愈的比例，称为恢复率系数。

在以上假设条件下，同时考虑脉冲接种因素，则对应的传染病动力学模型为：
/()()()/()()()()/()()()
/()()dS dt I t S t I t dE dt S t I t S t I t dI dt S t I t I t dR dt I t I t γτβββττβττδδγτ=--⎧⎪=---⎪⎨=---⎪⎪=--⎩（1）
()(1)() ()()()S nT p S nT n Z R nT R nT pS nT +-++--⎧=-∈⎨=+⎩（2）其中p 是类易感群体()S t 的脉冲接种率，T 为脉冲接种周期。

上述模型实质上是一个脉冲作用下具有时滞的微分方程组，对上述脉冲时滞微分模型进行数值仿真，就是对系统（1）
（2）求解，通过研究该方程组解的变化，从而得到传染病的发展趋势等相关内容。

3传染病动力学模型仿真系统的设计与实现
开发传染病动力学模型仿真系统的主要目的是建成一套能适应目前传染病动力学研究需要，且方便、快捷的数值仿真平台。

3．1系统组成
传染病动力学模型仿真系统主要分为四个部分：
1）模型分类系统可仿真的传染病动力学模型包括：常微分传染病模型、时滞微分传染病模型、偏微分传染病模型、常微分与偏微分混合型传染病模型等。

2）参数设定对模型中的各项参数进行设定，其中包括：对种群类别的设定（如仿真SIR 模型，即需选定易感类群体S (t )、染病类群体I (t )、恢复类群体R (t )）；仿真图形中曲线颜色、曲线线型以及曲线宽度、群体初始量的设定等。

此外还可以对传染病模型的相关系数进行设定，如：种群出生率、传染率系数、脉冲接种率、时滞量和垂直传染率等。

3）仿真图形显示系统图形仿真可将模型解的变化（即传染病的发展趋势等内容）以图像的形式显示出来，图像形式包括：二维曲线图，三维曲线图和三维曲面图。

4）文件输出系统可以将绘制的图形和数值试验数据以文件形式保存输出。

3．2系统采用的数值算法与仿真实现
考虑运算速度和精度的需要，系统对不同的模型采用不同的数值方法，其中常微分模型采用嵌入式Runge-Kutta 算法进行仿真，偏微分模型系统利用Matlab 中PDE 工具箱进行仿真。

下面只给出脉冲作用下非线性时滞传染病模型的数值解法。

首先时滞微分方程的一般形式如下：
00()(,(),()),()(),
y t f t y t y t t t y t t t t τφ'=-≥⎧⎨=≤⎩（3）这里：:[0,]N N N f C C C ∞⨯⨯→，:N Y R C →，0τ>，()N t C φ∈是给定向量函数，其中()y t τ-是方程的时滞项。

下面用Runge-Kutta 法给出求解上述问题的数值解法形式。

令(,,)A b c 表示给定的Runge-Kutta 方法，其中()ij s s A a ⨯=，12(,,,)T s b b b b = ，12(,,,)T
s c c c c = ，11s
i i b ==∑，[0,1]i c =，1i s ≤≤。

对应的Runge-Kutta 方法具体形式如下：
111111(,),1,,(,)s i n ij n j j j s n n i n j i i Y y h f t c h Y i s y y h b f t c h Y α--=--=⎧=++=⎪⎪⎨⎪=++⎪⎩∑∑ （4）
取步长/h m τ=，m 为正整数，将（4）式应用于方程（3），得：
()11()11(,,),1,,(,,)s n i n ij n j j j j s n n n i n j i j i Y y h f t c h Y Z i s y y h b f t c h Y Z α-=-=⎧=++=⎪⎪⎨⎪=++⎪⎩∑∑ （5）
如果用“~”符号表示精确解和数值解的近似，则~()n n y y t ，()~()n j n j Z y t c h τ+-代表第n 次迭代后时滞项的近似值，其中0n t t nh =+。

此外，对于脉冲项的处理，可以将离散点n t 整除脉冲周期T ，若n t 能整除脉冲周期T ，则n y 可加上或减去相应的数值。

将以上算法应用于系统（1）（2），其对应的MATLAB 核心代码如下：
%ST 等存贮离散网格点处计算值；tau 为时滞项；h 为步长%
ST=[];ET=[];IT=[];RT=[];m=tau/h;…
for i=1:tf/h %tf 为运算终值%
if (i<=m)%t 取(,0)τ-时，所取函数值%
z1=exp(t-tau);
z2=exp(2*(t-tau));
else %t >0时，所取离散网格点处数值%
z1=ST(1,i-m);
z2=IT(1,i-m);
end
％脉冲项处理，T 为脉冲接种周期％
if mod(i*h,T)==0
st=(1-p)*st;
rt=rt+p*st;
end
％以下为四级四阶Runge-Kutta 方法，F1，F2，F3为方程右端函数％
k11=F1(beta,gamma,y1,y2,z2);
k21=F2(beta,gamma,y1,y2,z1,z2);
k31=F3(beta,delta,y2,z1,z2);
k12=F1(beta,gamma,y1+h*k11/2,y2+h*k21/2,z2);
k22=F2(beta,gamma,y1+h*k11/2,y2+h*k21/2,z1,z2);
k32=F3(beta,delta,y2+h*k21/2,z1,z2);
k13=F1(beta,gamma,y1+h*k12/2,y2+h*k22/2,z2);
k23=F2(beta,gamma,y1+h*k12/2,y2+h*k22/2,z1,z2);
k33=F3(beta,delta,y2+h*k22/2,z1,z2);
k14=F1(beta,gamma,y1+h*k13,y2+h*k23,z2);
k24=F2(beta,gamma,y1+h*k13,y2+h*k23,z1,z2);
k34=F3(beta,delta,y2+h*k23,z1,z2);
st=st+h/6*(k11+2*k12+2*k13+k14);
et=et+h/6*(k21+2*k22+2*k23+k24);
it=it+h/6*(k31+2*k32+2*k33+k34);
rt=1-st-et-it;
ST=[ST st];ET=[ET y2];IT=[IT y3];RT=[RT y4];
t=t+h;
end
3．3仿真平台运行实例
以脉冲作用下的常微分SIR传染病模型为例，介绍传染病动力学模型仿真系统的部分界面。

图2为参数设定窗口，用户可以在该界面中选择模型的各项参数，包括群体类别、图形参数以及脉冲接种率等，图3为模型仿真结果的三维图形输出窗口。

图2选择模型参数窗口
图3模型仿真图像窗口
4结论
由于要建立与实际情况接近的数学模型，就需要增加模型系统的复杂性，从而对用纯粹数学方法研究模型造成了很大的困难。

如：种群规模变动的具有时滞的传染病模型的动力学性态、脉冲作用下具有年龄结构的AIDS模型等问题，目前在理论上还没有得以完全解决。

本文所开发的传染病动力学模型仿真系统以MATLAB为软件平台，应用GUI开发用户界面，并利用微分方程数值算法模拟传染病模型的各种性态，验证理论分析结果，从而为传染病动力学的相关科研工作提供了简捷、实用的仿真平台。

本文作者创新点：
本文利用MATLAB提供的图形化用户界面（GUI），结合时滞微分方程、脉冲微分方程等数值算法，开发了传染病动力学模型数值仿真平台，有效地解决了ODE45、DDE23等数值算法程序包无法直接对时滞微分传染病模型、脉冲传染病等新模型求解的问题。

参考文献
[1]M.E.Alexander,S.M.Moghadas.Periodicity in an epidemic model with a generalized nonlinear incidence[J].Math.Biosci.,2004,189(1):75–96.
[2]K.L.Cooke,P van den Driessche.Analysis of an SEIRS epidemic model with two delays.J. Math.Biol.,1996,35(2):240-260.
[3]魏巍,舒云星.基于非线性观测器的疾病发生率的估计[J].微计算机信息,2007,1:26-27.
[4]匡蛟勋.泛函微分方程的数值处理[M].北京:科学出版社.1999.。