基于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 引言

近年来，作为传染病研究的手段之一，利用计算机对传染病动力学模型进行数值仿真越来越受到人们的重视。诸如MA TLAB 中ODE45、DDE23等程序包，被人们普遍使用于传染病动力学模型的仿真中。近年来随着研究工作的深入，大量新的模型也逐渐受到人们的重视，如：时滞微分传染病模型；脉冲传染病模型；常微分、偏微分混合的传染病模型等。由于ODE45、DDE23等程序包不是针对传染病动力学模型所开发，无法解决以上这些模型的仿真问题，这些都给相关研究工作造成了一定的困难。本文利用MATLAB 提供的图形化用户界面（GUI ），结合时滞微分方程、脉冲微分方程等数值算法，并考虑传染病动力学模型的实际研究情况，开发了一套简单、实用的传染病动力学模型数值仿真平台。 2 传染病动力学模型的建立

从模型的数学结构来看，传染病动力学模型分为常微分模型、时滞微分模型、脉冲微分模型和偏微分模型等多种形式。以下以脉冲接种作用下的时滞传染病动力学模型为例，介绍模型的建立方法。“时滞”可以反映传染病的潜伏期，患者对疾病的感染期和恢复者对疾病的免疫期等实际现象，因此使用“时滞”模型更贴近实际。如Cooke 等人将时滞因素引入到SEIRS 传染病模型中，用时滞项来反映传染病的潜伏期，建立了如图1所示的仓室框图。

图1 SEIRS 模型的仓室框图

在此模型中，将传染地区的人群分为四类：用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 = ，1

1s

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）