一類復雜流行病學模型的混沌研究*

狄根虎 許 勇 徐 偉 顧仁財

(西北工業大學應用數學系，西安 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.