一类复杂流行病学模型的混沌研究*

狄根虎 许 勇 徐 伟 顾仁财

(西北工业大学应用数学系，西安 710129)

一类复杂流行病学模型的混沌研究*

狄根虎 许 勇徐 伟 顾仁财

(西北工业大学应用数学系，西安 710129)

(2010年5月15日收到;2010年9月17日收到修改稿)

研究了一类周期变化的非线性复杂发病率的广义流行病学模型(SIR(susceptible，infected，recovered)模型).通过一系列坐标变换将原模型转化为Hamilton系统，运用Melnikov方法证明了该系统存在混沌运动，给出了发生同宿分岔的条件，并用数值模拟验证了上述结果.

SIR(susceptible，infected，recovered)模型，混沌运动，Melnikov方法，同宿分岔

PACS:05.45.-a，05.45.Pq，02.30.Qz

1.引 言

数学模型和方法已经用于传染病的传播原理研究.Kermack和 McKendrick［1］在研究1665—1666年黑死病的流行规律时，构造了著名的((SIR)仓室模型，现今最常用的是Anderson和May［2］所建立的传染病数学模型，即假设个体对传染病是免疫情形下的标准SIR模型:

其中，模型(1)假定人口是不变的，将人口分为三类:S(t)代表t时刻易感人群;I(t)代表t时刻的感染人群;R(t)代表t时刻的恢复且不再感染这类疾病的人群.参数μ表示人口出生率或死亡率(此模型假设出生率和死亡率相等)，γ表示感染人群的恢复率，而B(I，S，t)代表随时间周期变化的发病率函数.

在以往研究中，通常取B(I，S，t)=β(t)IS，其中传染率β(t)或者是常数或者周期变化，如文献［3—5］中，取β(t)=β0(1+β1sin(ωt))，并数值证明了当 β1充分大时，该模型存在混沌运动.文献［6—10］研究了复杂发病率函数，得到了具有不同于双线性发病率函数的动力学行为，包括霍普夫分岔、鞍结分岔和同宿分岔.这些分岔意味着传染病的爆发和灭绝，由参数 β，p，q决定.文献［11—17］基于复杂网络对一类传染病模型进行研究.因此，具有复杂发病率函数的SIR模型(1)的研究是重要和必要的.然而，具有复杂发病率函数的 SIR模型(1)的混沌行为还研究较少，本文取复杂发病率函数为B(I，S，t)=β(t)IpSq，传染率函数取为

其中，b1，b2，b3，Ω以及 p，q都为常数且 p>1，使得2pb2-(p-1)b21>0，ε>0是小参数.这里我们所取的传染率函数在 p=2，q=1时就退化到文献［18］中所考虑的传染率情形，它是在文中第二部分的模型变换时，为了将模型变换为Hamilton系统的标准形式，对传染率函数进行了适当的修正，考虑了参数p，q对传染率的影响，这也是合理的.同时假定模型(1)中的恢复率函数为γ=μ(1+b1ε2)/(p-1).

本文将在复杂发病率函数(2)下研究模型(1)的混沌动力学行为.首先通过一系列坐标变换将模型(1)化为Hamilton系统，然后应用Melnikov方法计算了混沌产生的临界值，得到同宿分岔的条件，并通过数值实验说明了解析结果的正确性.

2.模型变换

假定模型(1)中的总人口是常数，且可表示为S +R+I=1.

将S=1-I-R代入模型(1)，模型简化为

由于μ>0，定义新时间变量 t′=μt，消去 μ，方程(3)为

其中，β=μβ′，γ=μγ′.

令y=γ′I-R，即 I=γ′(-1)(y+R)，代入(4)式有

其中，(5)式中的导数是对新的时间变量t′而言.现在计算(5)式稳态点.当y=0时，有

显然R=0是方程(6)的一个解，即具有稳态点(0，0)，此时称(S，I，R)=(1，0，0)为无病平衡点.进一步，方程(6)简化为

记 H=γ′/(1+γ′)，σ=(1+γ′)/β′γ′(1-p)，则Rp-1(1-R/H)q=σ，R∈［0，H］

定义函数f(R)=Rp-1(1-R/H)q，其中，R∈［0，H］.由极值原理，可知f(R)在［0，H］上的最大值为 σ*，其中，σ*=(p-1)p-1qqHp-1/(p+q-1)p+q-1.

当σ>σ*时，系统(7)的正平衡点不存在，当σ=σ*时，存在一个正平衡点，当 σ<σ*时，存在两个正平衡点.这里我们给出如下定理:

定理 如果β′<(p+q-1)p+q-1(1+γ′)p/((p-1)p-1qq)，系统(7)除无病平衡点外不存在正平衡点;如果 β′=(p+q-1)p+q-1(1+γ′)p/ ((p-1) q)，系统将除无病平衡点外还存在一个正平衡点;如果β′>(p+q-1)p+q-1(1+γ′)p/((p-1)p-1qq)，系统除无病平衡点外还存在两个平衡点.

由上述定理可得，系统(7)在β′=βs′n=(p+q-1)p+q-1(1+γ′)p/((p-1)p-1qq)处发生鞍结分岔，在此处的平衡点为

系统(7)在稳态(R，0)处的Jacobi矩阵为

其中，

注意到发生鞍结分岔时，J21=0.如果J22=0，则稳态点就处在更加退化的情形，此时将稳态点称为Takens-Bogdanov点，由［19］知，在 Takens-Bogdanov点附近可以将方程化为近似 Hamilton系统的标准形式.

在鞍结分岔时，如果γ′=1/(p-1)，则矩阵A的迹消失，即J22=0.

取上述参数值，可得退化平衡点为

由于霍普夫分岔的曲线是在此退化平衡点终止的，曲线是由 J21<0，J22=0给定的.为扰动退化点，定义新变量和新参数

其中，x，g，b表示正的小扰动量.

将其带入(5)式，忽略x，g，b的三次项及高阶项，有

其中，

引进小参数ε>0(与(2)式中的小参数ε相同)，令

将(13)式代入(12)式，忽略高次项得

定义新的坐标变量和时间变量x=ε2u，y=ε3v，τ=εt′，系统(11)为

令z=u-((p-1)2/(p2(p+q-1)))b1，则系统(15)化为

定义变量k=(p(p+q-1)/(2q(1-p)))z，w=(p(p+q-1)/(2q(1-p)))v，方程(16)为

其中，c2=(p2(p+q-1)/(2(p-1)q))b2- (p(p+q-1)/(4q)，由前面假设可知 c2是严格大于0的常数.以后取c>0.

当b3=0，方程(17)是Takens-Bogdanov分岔的标准形式.假定 b1>0，b3=0以及ε是充分小的正常数，则方程有两平衡点S±=(±c，0)，S+是一鞍点，S-有如下三种情形:

3.SIR模型中的混沌

本节运用 Melnikov方法［20—24］研究 SIR模型经过变换后的Hamilton系统(17)，证明系统的混沌存在性以及同宿分岔的条件.

如果ε=0，系统(17)的未扰形式为

其中，导数是对新时间变量τ而言，且c>0.

系统(18)的Hamilton函数和势函数分别为

方程(18)有两平衡点，(c，0)是鞍点，(-c，0)是中心.中心是被连续的一族周期轨道(kT(τ)，wT(τ))所围绕，而且这些周期轨道是有界限的，是被同宿轨道(kh(τ)，wh(τ))所控制，这是因为

系统的Melnikov函数为

未扰系统的同宿解为［25］

代入(19)式并积分得

令b1=0，给出一简单表达形式，其中

令b3=0，由(21)式给出系统无周期激励的同宿分岔条件为b1=10(p-1)c/(7p)，即同宿分岔是由参数空间(b1，b2)上的曲线

确定的.

4.数值结果

本节运用数值实验来验证 Melnikov方法所得到的解析结果，即证实系统存在混沌运动.针对系统(17)，取如下参数:

图1给出了上述参数情形下的Poincaré截面，而图1中点的杂乱分布，由 Poincaré理论知系统具有混沌运动，与解析结果是一致的.图2给出了系统(17)的相图，我们发现系统也是混沌的.很明显，数值实验与解析结果是一致的.

图1 系统(17)的Poincaré截面 (a)bM=2.7919，(b)bM=6.6050

图2 系统(17)的相图 (a)bM=2.7919，(b)bM=6.6050

5.结 论

本文通过解析方法即 Melnikov方法证明复杂周期激励的传染病模型在 Takens-Bogdanov分岔点的很小范围内具有混沌运动，给出出现混沌的临界表达式(22)式.数值模拟验证了解析结果的正确性.从图1，2可知，参数达到临界值时，传染病会处于爆发状态，因此在预防传染病的传播时，要对参数进行合理控制，使得传染病在可控范围内.

［1］Kermack W O，McKendrick A G 1927 Proc.Roy.Soc.A 115 700

［2］Anderson R M，May R M 1991 Infectious Diseases of Humans (Oxford:Oxford University Press)

［3］Aron J L，Schwartz I B 1984 Theor.Biol.110 665

［4］Kuznetsov，Yu A，Piccardi C 1994 Math.Biol.32 109

［5］Olsen L F，Shaffer W M 1990 Science 249 499

［6］Hethcote H W，Driessche P V D 1991 Math.Biol.29 271

［7］Liu W M，Hethcote H W，Levin S A 1987 Math.Biol.25 359

［8］Liu W M，Levin S A，Iwasa Y 1986 Math.Biol.23 187

［9］Yi H T 2009 MS dissertation(Hunan:Hunan university)(in Chinese)［尹海涛2009硕士学位论文(湖南:湖南大学)］

［10］Li J 2005 MS dissertation(Nanjing:Nanjing university of science and technology)(in Chinese)［李 静 2005硕士学位论文(南京:南京理工大学)］

［11］Ni S J，Weng W G，Fan W C 2009 Acta Phys.Sin.58 3707(in Chinese)［倪顺江、翁文国、范维澄2009物理学报58 3707］

［12］Liu M X，Ruan J 2009 Chin.Phys.B 18 2115

［13］Liu M X，Ruan J 2009 Chin.Phys.B 18 5111

［14］Shi H J，Duan Z S，Chen G R，Li R 2009 Chin.Phys.B 18 3309

［15］Pei W D，Chen Z Q，Yuan Z Z 2008 Chin.Phys.B 17 0373

［16］Zhong L，Weng J Q 2000 Acta Phys.Sin.49 0626(in Chinese)［钟 玲、翁甲强2000物理学报49 0626］

［17］Xu D，Li X，Wang X F 2007 Acta Phys.Sin.56 1313(in Chinese)［许 丹、李 翔、汪小帆2007物理学报56 1313］

［18］Glendinning P，Perry L P 1997 Math.Biol.34 359

［19］Chow S N，Li C Z，Wang D 1994 Normal forms and bifurcation of plannarvectorfields(Cambridge:Cambridge University Press)

［20］Yang X L，Xu W，Sun Z K 2006 Acta Phys.Sin.55 1678(in Chinese)［杨晓丽、徐 伟、孙中奎2006物理学报55 1678］

［21］Liu Z R 2002 The analytical method of chaos(Shanghai: Shanghai university press)p56(in Chinese)［刘增荣 2002混沌研究中的解析方法(上海:上海大学出版社)第56页］

［22］Lei Y M，Xu W 2007 Acta Phys.Sin.56 5103(in Chinese)［雷佑銘、徐 伟2007物理学报56 5103］

［23］GuckenheimerJ， Holmes P 1983 NonlinearOscillations，Dynamical Systems and Bifurcations of Vectors Fields(New York: New York Springer)

［24］Yang X L，Xu W 2009 Acta Phys Sin.58 3722(in Chinese)［杨晓丽、徐 伟2009物理学报58 3722］

［25］Baesens N C，Nicolis G 1983 Z.Phys.B 52 345

PACS:05.45.-a，05.45.Pq，02.30.Qz

Chaos for a class of complex epidemiological models*

Di Gen-Hu Xu YongXu Wei Gu Ren-Cai
(Department of Applied Mathematics，Northwestern Polytechnical University，Xi'an 710129，China)

15 May 2010;revised manuscript

17 September 2010)

We study the well-known SIR(susceptible，infected，recoverd)model with nonlinear complex incidence rates. Firstly，a series of coordinate transformations are carried out to change the equations as the amenable Hamiltonian systems. Secondly the Melnikov′s method is used to establish the conditions of existence of chaotic motion and find the analytically critical values of homoclinic bifurcation.Good agreement can be found between numerical results and analytical results.

SIR(susceptible，infected，recoverd)model，chaotic motion，Melnikov's method，homoclinic bifurcation

*国家自然科学基金(批准号:10972181，11002001，10872165，10932009)、西北工业大学基础研究基金和翱翔之星资助的课题.

*Project supported by the National Natural Science Foundation of China(Grant Nos.10972181，11002001，10872165，10932009)，the NPU Foundation for Fundamental Research and the“Aoxiang Star”Plant of Northwestern Polytechnical University，China.