一類具有時(shí)滯和進(jìn)化效應(yīng)的SIR模型的穩(wěn)定性和Hopf分支分析

摘 要：本文主要研究了一類在齊次Neumann邊界條件下具有時(shí)滯和進(jìn)化效應(yīng)的SIR模型.首先，以時(shí)滯為參數(shù)，分析了時(shí)滯對(duì)該系統(tǒng)正平衡點(diǎn)穩(wěn)定性的影響，并給出了Hopf分支的存在條件.其次，利用規(guī)范型理論和中心流形定理，得到Hopf分支方向和分支周期解的穩(wěn)定性.最后，利用Matlab進(jìn)行數(shù)值模擬，驗(yàn)證結(jié)論的正確性，得出時(shí)滯會(huì)使穩(wěn)定的系統(tǒng)變得不穩(wěn)定，并產(chǎn)生Hopf分支.

關(guān)鍵詞：時(shí)滯； Hopf分支； 穩(wěn)定性； SIR模型

中圖分類號(hào)：O175.26

文獻(xiàn)標(biāo)志碼： A

Stability and Hopf bifurcation analysis of the SIR model with time-delay and evolutionary effects

YUAN Hai-long， WANG Jia-yue

（School of Mathematics and Data Sciencnce， Shaanxi University of Science amp; Technology， Xi′an 710021， China

）

Abstract：In this paper，we study a class of SIR model with time-delay and evolutionary effects under the homogeneous Neumann boundary conditions.Firstly，the stability of the positive equilibrium point is analyzed with time delay as the parameter，and the existence condition of Hopf bifurcation is given.Secondly，the bifurcation direction and the stability of periodic solutions are given by using the theory of normal form and center manifold.Finally，Matlab was used for the numerical simulation to verify the correctness of the conclusion，and it is concluded that the time delay will make the stable system become unstable and generate Hopf bifurcation.

Key words：time-delay；" Hopf bifurcation； stability； SIR mode

0 引言

自1927年Kermack與Mckendrick給出了傳染病倉(cāng)室模型［1，2］，建立了經(jīng)典的閾值理論，微分方程和動(dòng)力系統(tǒng)理論開始被廣泛用于研究多種疾病的起源、進(jìn)化以及傳播.傳染病動(dòng)力學(xué)模型由此也得到了迅速的發(fā)展，這些模型在不同的傳染病背景下給出了疾病傳播的一般性規(guī)律，使人們對(duì)相關(guān)疾病的傳播可以有更深刻的理解，其結(jié)果有助于預(yù)測(cè)傳染病的發(fā)展趨勢(shì)，確定疾病傳播的關(guān)鍵因素，尋找預(yù)防和控制傳染病傳播的最佳策略.在此研究過程當(dāng)中，許多學(xué)者都利用SIR模型研究傳染病的傳播規(guī)律［3-6］.

本文研究以健康個(gè)體、已感染個(gè)體和已恢復(fù)個(gè)體為研究對(duì)象的SIR模型，是基于文獻(xiàn)［7］中研究的模型，進(jìn)一步考慮了在時(shí)滯和空間異質(zhì)以及疾病的傳染力是具有飽和狀態(tài)的基礎(chǔ)上，由于健康個(gè)體自身的防御成本和已感染個(gè)體感染能力有適應(yīng)性變化，從而對(duì)系統(tǒng)動(dòng)力學(xué)行為產(chǎn)生進(jìn)化效應(yīng)影響的情況［8-10］，所以此問題更加貼近生活，更具現(xiàn)實(shí)意義.

由于健康個(gè)體的內(nèi)稟增長(zhǎng)率取決于自身對(duì)病毒的防御成本，還有感染率由健康人群防御成本和已感染者感染能力這兩大因素，文獻(xiàn)［11］考慮了具有進(jìn)化效應(yīng)的模型：

dSdt=r（u）（1-S+IK）S-β（u，v）SIS+I，

dIdt=β（u，v）SIS+I-（δ+γ）I，

dRdt=γI-μR

（1）

而在實(shí)際生活中，一個(gè)患者與他人的接觸能力總歸是有限的，所以疾病的傳染力有一個(gè)飽和狀態(tài).因此，研究飽和發(fā)生率的傳染病模型也具有一定的現(xiàn)實(shí)意義.Anderson等［12］提出了用飽和發(fā)生率的數(shù)學(xué)模型來(lái)分析流行病的性態(tài)，Ruan等［3］引入形如βSIP1+αIP的飽和傳染率，并分析了模型的分支情況.在考慮P=1的情況下，系統(tǒng)（1）可被改進(jìn)為如下模型：

dSdt=r（u）（1-S+IK）S-β（u，v）SI1+αI，

dIdt=β（u，v）SI1+αI-（δ+γ）I，

dRdt=γI-μR

（2）

由于所有傳染病都有潛伏期，因此將感染媒介在病媒中發(fā)展的潛伏期時(shí)間引入流行病模型［6，13］會(huì)更符合實(shí)際.比如在潛伏期的個(gè)體雖已被感染，但沒有傳染性［13］，病原生物體在被感染的宿主體內(nèi)繁殖需要足夠多的時(shí)間以后，才對(duì)其他宿主具有傳染性［6］，因此具有時(shí)滯的流行病模型逐漸引起了廣泛的關(guān)注［14］.在許多種群動(dòng)力學(xué)模型中，時(shí)滯可以使平衡變得不穩(wěn)定，從而導(dǎo)致Hopf分支的產(chǎn)生［15］.鑒于這一點(diǎn)，假設(shè)潛伏期τ是恒定的，并假設(shè)t時(shí)刻的感染力依賴于前一個(gè)t-τ時(shí)刻的易感人群和具有傳染性的人群，則具時(shí)滯的流行病模型如下：

dSdt=r（u）（1-S+IK）S-β（u，v）S（t-τ）I（t-τ）1+αI（t-τ），

dIdt=β（u，v）S（t-τ）I（t-τ）1+αI（t-τ）-（δ+γ）I，

dRdt=γI-μR

（3）

此外，在構(gòu)建現(xiàn)實(shí)的流行病模型時(shí)，隨著人們對(duì)空間方面的認(rèn)識(shí)越來(lái)越多［16-18］，人們發(fā)現(xiàn)易感人群和具有傳染性的人群會(huì)在空間域上隨機(jī)移動(dòng).若要將該運(yùn)動(dòng)過程考慮進(jìn)流行病模型，其常見方法是假設(shè)某些宿主是隨機(jī)運(yùn)動(dòng)的，并據(jù)此建立反應(yīng)擴(kuò)散方程［19］.假設(shè)易感人群和感染人群是隨機(jī)移動(dòng)的，那么在有界開區(qū)間（0，lπ），l∈R+上考慮模型：

St=dSΔS+r（u）（1-S+IK）S-β（u，v）SτIτ1+αIτ，

x∈（0，lπ），tgt;0，

It=dIΔI+β（u，v）SτIτ1+αIτ-（δ+γ）I，x∈（0，lπ），tgt;0，

Rt=dRΔR+γI-μR，x∈（0，lπ），tgt;0，

νS=νI=νR=0，x=0，lπ，tgt;0，

S（x，0）=S0（x）≥0，I（x，0）=I0（x）≥0，

R（x，0）=R0（x）≥0，x∈＼，t∈-τ，0

（4）

其中，S（x，t）、I（x，t）和R（x，t）代表易感人群、感染人群和恢復(fù)人群在位置x和時(shí)間t處的密度，正常數(shù)dS、dI、dR分別為易感人群、感染人群和恢復(fù)人群的擴(kuò)散系數(shù)，δ為已感染個(gè)體的死亡率，μ為已恢復(fù)個(gè)體的死亡率，γ為個(gè)體的恢復(fù)率，r0為最大增長(zhǎng)率，c1為防御成本，K為環(huán)境最大容納量，β0為最大感染率，θ為感染率在不同人群之間差異的靈敏度，為了滿足實(shí)際生物意義，S（0）+I（0）lt;K.另外，關(guān)于r（u），β（u，v）的假設(shè)如下［11］：

（1）健康個(gè)體防御成本u和健康個(gè)體的內(nèi)稟增長(zhǎng)率r之間的關(guān)系假設(shè)為r（u）=r0e-c1u2（參見文獻(xiàn)［8］），因此隨著人群總數(shù)接近環(huán)境最大容納量，r對(duì)動(dòng)力學(xué)的影響減小.當(dāng)健康人群占有率低時(shí)，相同水平的防御成本較高，而當(dāng)健康人群占有率高時(shí)，相同水平的防御成本較低.

（2）易感人群的感染率為β（u，v）=β01+eθ（u-v），

下文將用β代替β（u，v），其中v為已感染個(gè)體的傳染能力，這是由防御成本和感染能力決定的［9］，感染率隨著防御成本的增加而減小，隨著感染能力的增強(qiáng)而增強(qiáng).

本文研究了時(shí)滯在偏微分系統(tǒng)中對(duì)正平衡點(diǎn)產(chǎn)生的影響，并分析了該平衡點(diǎn)的穩(wěn)定性和Hopf分支的存在性，給出了判斷分支方向和分支周期解穩(wěn)定性的表達(dá)式.從分析中可以得出，當(dāng)時(shí)滯參數(shù)小于某一臨界值時(shí)，正平衡點(diǎn)的穩(wěn)定性不發(fā)生改變.當(dāng)時(shí)滯參數(shù)大于某一臨界值時(shí)，會(huì)改變正平衡點(diǎn)的穩(wěn)定性，系統(tǒng)在其附近產(chǎn)生振蕩并產(chǎn)生Hopf分支.

本文主要工作內(nèi)容分為三個(gè)部分.第1部分以時(shí)滯為分支參數(shù)，通過研究正平衡點(diǎn)對(duì)應(yīng)特征方程，得到平衡點(diǎn)的穩(wěn)定性結(jié)果以及Hopf分支的存在性；第2部分運(yùn)用中心流形定理和規(guī)范型理論得到Hopf分支的分支方向和分支周期解的穩(wěn)定性； 第3部分利用Matlab軟件進(jìn)行數(shù)值模擬來(lái)驗(yàn)證本文的理論分析結(jié)果.本文用N0表示非負(fù)整數(shù)集，R+表示所有正實(shí)數(shù)的集合.

1 穩(wěn)定性與分支分析

由于系統(tǒng)（1）中的第三個(gè)方程與前兩個(gè)方程解耦，因此系統(tǒng)可被簡(jiǎn)化為如下系統(tǒng)：

St=dSΔS+r（u）（1-S+IK）S-βS（t-τ）I（t-τ）1+αI（t-τ），

x∈（0，lπ），tgt;0，

It=dIΔI+βS（t-τ）I（t-τ）1+αI（t-τ）-（δ+γ）I，x∈（0，lπ），tgt;0，

νS=νI=νR=0，x=0，lπ，tgt;0，

S（x，0）=S0（x）≥0，I（x，0）=I0（x）≥0，

R（x，0）=R0（x）≥0，x∈＼，t∈［-τ，0］

（5）

通過文獻(xiàn)［11］可知，

（1）當(dāng)R0lt;1時(shí)，系統(tǒng)（5）存在無(wú)病平衡點(diǎn)E1（K，0）;

（2）當(dāng)R0gt;1，（δ+γ）lt;βlt;r（u）+（δ+γ）時(shí)，系統(tǒng)（5）存在唯一正平衡點(diǎn)E*（S*，I*），且當(dāng)給u，v賦不同值時(shí)，β，r（u）為常數(shù)，其中

S*=（r（u）-β+（δ+γ））（δ+γ）Kr（u）β，

I*=S*（β-（δ+γ））δ+γ.

下面對(duì)正平衡點(diǎn)E*（S*，I*）做變量替換，則系統(tǒng)（5）變?yōu)椋?/p>

tS=dSΔS+r（u）（S+S*）-r（u）K（S+S*）2

-r（u）K（S+S*）（I+I*）

-β（Sτ+S*）（Iτ+I*）1+α（Iτ+I*），x∈（0，lπ），tgt;0，

tI=dIΔI+β（Sτ+S*）（Iτ+I*）1+α（Iτ+I*）

-（δ+γ）（I+I*），x∈（0，lπ），tgt;0，

νS=νI=νR=0，x=0，lπ，tgt;0，

S（x，0）=S0（x）≥0，I（x，0）=I0（x）≥0，

R（x，0）=R0（x）≥0，x∈＼，t∈［-τ，0］

（6）

其中，S=S（x，t），I=I（x，t），Sτ=S（x，t-τ），Iτ=I（x，t-τ），那么E*（S*，I*）的穩(wěn)定性問題就轉(zhuǎn)換為原點(diǎn)（0，0）的穩(wěn)定性問題.

定義X=C（［0，lπ］，R2），在空間C（［-τ，0］，X］中，系統(tǒng)（6）有如下所示的抽象泛函微分方程表達(dá)式：

dU（t）dt=DΔU（t）+L（Ut）+F（Ut）

（7）

其中，DΔ=（d1Δ，d2Δ），dom（DΔ）=（S，I）T：S，

I∈C2（0，lπ，R），Sx，Ix=0，x=0，lπ，且L：

C（-τ，0，X）→X，F(xiàn)：C（-τ，0，X）→X，對(duì)

于=（1，2）T∈C（-τ，0）有

L（）=r（u）1-2S*+I*K1（0）-βI*1+αI*1（-τ）

-r（u）KS*2（0）-βS*（1+αI*）22（-τ）

βI*1+αI*1（-τ）-（δ+γ）2（0）+βS*（1+αI*）22（-τ），

F（）=F1（）

F2（），

其中

F1（）=r（u）S*-r（u）K（12（0）+S*2+

1（0）2（0）+S*I*）-

β（1（-τ）+S*）（2（-τ）+I*）1+α（2（-τ）+I*）+

β（1+αI*）1（-τ）I*+βS*2（-τ）（1+αI*）2，

F2（）=β（1（-τ）+S*）（2（-τ）+I*）1+α（2（-τ）+I*）-

β（1+αI*）1（-τ）I*+βS*2（-τ）（1+αI*）2-

（δ+γ）I*.

則系統(tǒng)（6）在（0，0）處的線性化系統(tǒng)為

dU（t）dt=DΔU（t）+L（Ut）

（8）

其中，U（t）=（1（t）

2（t）），D=（d10

0d2）.

從Wu［20］中，可以得到式（8）的特征方程為

-λy+DΔy+L（eλ·y）=0，y∈dom（DΔ），y≠0

（9）

由特征值問題

-Δφ=λφ，x∈（0，lπ），φ′（0）=φ′（lπ）=0，

可得特征值λn=n2l2（n∈N0），其相應(yīng)的特征函數(shù)為φn（x）=cosnlx. 將y=∑∞n=0cosnlxy1n

y2n代入式（9）得到

-d1n2l2+m-we-λτ-p-qe-λτ

we-λτ-d2n2l2-z+qe-λτy1n

y2n=

λy1n

y2n，n∈N0

其中m=r（u）1-2S*+I*K，w=βI*1+αI*，

p=r（u）KS*，q=βS*（1+αI*）2，z=（δ+γ），

則式（9）等價(jià)于

Δn（λ，τ）=λ2+Anλ+Bnλe-λτ+Cne-λτ+Dn=0，n∈N0

（10）

其中

An=d1n2l2-m+d2n2l2+z，Bn=w-q，

Cn=-d1n2l2q+mq+d2n2l2w+wz+pw，

Dn=d1d2n4l4+zd1n2l2-md2n2l2-mz.

由于E*（S*，I*）的穩(wěn)定性是由式（10）根的分布來(lái)決定的，則當(dāng)式（10）所有的根具有負(fù)實(shí)部時(shí)，E*（S*，I*）局部漸近穩(wěn)定.

首先討論系統(tǒng)（5）中τ=0時(shí)，式（10）根的分布情況.

引理1 對(duì)于系統(tǒng)（5），當(dāng)τ=0，且滿足R0gt;1，（δ+γ）lt;βlt;r（u）+（δ+γ）時(shí)，n∈N0，

Δn（λ，0）=0的根都具有負(fù)實(shí)部，平衡點(diǎn)E*（S*，I*）是局部漸近穩(wěn)定的.

證明：當(dāng)τ=0時(shí)，式（10）變成

Δn（λ，0）=λ2+（An+Bn）λ+（Cn+Dn）=0，

n∈N0

（11）

由于n→∞時(shí)，

An+Bn=d1n2l2+d2n2l2-m+z+w-q→∞

Cn+Dn=d1d2n4l4+（d1z-d1q-d2m+d2w）n2l2+

m（q-z）+w（z+p）→∞.

所以，δlt;0，使得sup{λ：n∈N0，Δn（λ，0）=0}lt;δ， 則當(dāng)τ=0時(shí)，平衡點(diǎn)E*（S*，I*）是局部漸近穩(wěn)定的.

證畢.

接下來(lái)，討論系統(tǒng)（5）中τgt;0時(shí)，式（10）根的分布情況.

若λ=0是式（10）的根，則有

Cn+Dn=d1d2n4l4+（d1z-d1q-d2m+d2w）n2l2+

m（q-z）+w（z+p）=0，

由引理1知，Bn+Cgt;0，所以λ=0不是式（10）的根.

引理2 對(duì)于式（10），n∈N0，Δn（λ，τ）=0

沒有純虛根.

證明：假設(shè)±iω（ωgt;0）是式（10）的一對(duì)純虛根，將λ=iω（ωgt;0）代入有

-ω2+Aniω+Bniωe-iωτ+Cne-iωτ+Dn=0，

分離其實(shí)部與虛部可得

ω2-Dn=ωBnsinωτ+Cncosωτ，

Anω=Cnsinωτ-ωBncosωτn∈N0

（12）

有

ω4+ω2（An2-2Dn-Bn2）-Cn2+Dn2=0，n∈N0

（13）

其中

An2-2Dn-Bn2=n4l4（d12+d22）-n2l2（3d1m+d2m-2d1z-2d2z）+m2+z2-3mz-（w-q）2，

Dn2-Cn2=d12d22n8l8+n6l6（2d12d2z-2md1d22）

+n4l4（d1z-d2m）2-2mzd1d2

-（d1q-d2w）2+n2l2［2mz2（d2-d1）+

2d1（mq2+wzq+pwq）-2d2（w2z+pw2+wmq）］

+m2（z2+q2）+w2（z2+p2）+2mqw（z+p）.

令z=ω2，帶入式（13）得：

ω2+ω（An2-2Dn-Bn2）+Dn2-Cn2=0，n∈N0

（14）

且limn→∞（An2-2Dn-Bn2）gt;0，limn→∞（Dn2-Cn2）gt;0，則式（14）沒有正根，式（10）沒有純虛根.

證畢.

接下來(lái)研究式（14）存在正根的條件.對(duì)于式（11）有l(wèi)imn→0（D2n-C2n）lt;0，limn→∞（D2n-C2n）gt;0，則存在一個(gè)最小非負(fù)數(shù)N，使得式（13）在ngt;N的情況下無(wú)正根，在0≤n≤N的情況下最多有一個(gè)正根，且該正根ωn滿足

ωn2=-（An2-2Dn-Bn2）2

+（An2-2Dn-Bn2）2-4（Dn2-Cn2）2

（15）

且由式（12）可得

F（ω）=sinωτ=Bnωn3-BnDnωn+AnCnωnCn2+Bn2ωn2，

G（ω）=cosωτ=ωn2（Cn-AnBn）-CnDnCn2+Bn2ωn2，

則式（10）有一對(duì)純虛根，當(dāng)

τ=τjn=τ0n+2jπωn，j∈N0

（16）

其中τon滿足

τ0n=arccos（F（ωn））ωn，""""""" G（ωn）≥0，

2π-arccos（F（ωn））ωn，G（ωn）lt;0

（17）

引理3 令λn（τ）=γn（τ）+iωn（τ）為式（10）的根，當(dāng)τ趨近于τjn時(shí)，滿足γn（τjn）=0，ωn（τjn）=ωn，此時(shí)橫截條件得以滿足，即當(dāng)0≤n≤N，j∈N0，有γn′（τjn）gt;0成立.

證明：將λn（τ）帶入式（10）并對(duì)其兩邊同時(shí)對(duì)τ求導(dǎo)得

dλndττ=τjn-1=（2λn（τ）+An）eλn（τ）τjn+Bnλn（τ）（Bnλn（τ）+Cn）-τjnλn（τ），

又因?yàn)棣豱和τjn滿足式（12）和式（15），帶入上式得

Redλndττ=τjn-1=2Bnωn2sinωnτjn+2CnωncosωnτjnBn2ωn2+Cn2

+AnCnsinωnτjn-AnBnωncosωnτjn-Bn2ωnBn2ωn2+Cn2

=2Bnωn2sinωnτjn+2CnωncosωnτjnBn2ωn2+Cn2

+AnCnsinωnτjn-AnBnωncosωnτjn-Bn2ωnBn2ωn2+Cn2

=2ωn（Bnωnsinωnτjn+Cncosωnτjn）Bn2ωn2+Cn2

+An（Cnsinωnτjn-Bnωncosωnτjn）-Bn2ωnBn2ωn2+Cn2

=2ωn（ωn2-Dn）+An2ωn-Bn2ωnBn2ωn2+Cn2

=ωn2ωn2+（An2-2Dn-Bn2）Bn2ωn2+Cn2gt;0.

因此，橫截條件成立.

證畢.

顯然，根據(jù)式（16）可得τjn關(guān)于j是單調(diào)遞增的，即τj+1ngt;τjn.接下來(lái)討論τjn關(guān)于n的單調(diào)性.

引理4 對(duì)于R0gt;1，（δ+γ）lt;βlt;r（u）+（δ+γ）且0≤n≤N，j∈N0，有ω2n+1lt;ω2n成立.

證明：對(duì)式（15）作恒等變形得

ωn2=-（An2-2Dn-Bn2）+（An2-2Dn-Bn2）2-4（Dn2-Cn2）2=

2（An2-2Dn-Bn2）2（Cn2-Dn2）2+4Cn2-Dn2+An2-2Dn-Bn2Cn2-Dn2.

由于An2-2Dn-Bn2關(guān)于n嚴(yán)格遞增，且Cn2-Dn2關(guān)于n嚴(yán)格遞減，其中0≤n≤N，則有ωn+12lt;ωn2.

證畢.

下面給出τjn關(guān)于n的單調(diào)性.

引理5 對(duì)于R0gt;1，（δ+γ）lt;βlt;r（u）+（δ+γ）且0≤n≤N，j∈N0，有τjn+1gt;τjn.

證明：由式（15）得τ0n=arccosF（ωn）ωn， 進(jìn)一步可得τ0n+1gt;τ0n，其中（0≤n≤N）.由于ωn+1lt;ωn，則由式（14）得τjn+1gt;τjn，（j≥1，0≤n≤N）.

證畢.

根據(jù)引理6可知，系統(tǒng)（6）的第一個(gè)Hopf分支值為τ00.綜合上述分析，得到以下定理.

定理1 對(duì)于系統(tǒng)（6），

（1）如果τ∈［0，τ00），則式（10）的所有根都具有負(fù)實(shí)部，即地方病平衡點(diǎn)E*（S*，I*）是局部漸近穩(wěn)定的.

（2）如果τ=τjn（0≤n≤N，j∈N0），則式（10）除±iω之外的所有根都具有負(fù)實(shí)部，并且此時(shí)系統(tǒng)（6）在E*（S*，I*）處經(jīng)歷Hopf分支.

（3）如果τgt;τ00，則式（10）至少有兩個(gè)根具有正實(shí)部，地方病平衡點(diǎn)E*（S*，I*）變得不穩(wěn)定.

2" Hopf分支穩(wěn)定性及其分支方向

在第1節(jié)中，已得到當(dāng)τ=τjn時(shí)，系統(tǒng)（5）在E*（S*，I*）處產(chǎn)生Hopf分支.本節(jié)將運(yùn)用中心流形定理和規(guī)范型理論［20，21］研究當(dāng)τ=τ0=τ00時(shí)，Hopf分支的方向和分支周期解的穩(wěn)定性.

假設(shè)τ0為系統(tǒng)（5）的第一個(gè)Hopf分支值，令τ=τ0+μ，則當(dāng)μ=0時(shí)，τ是系統(tǒng)（6）的Hopf分支值，再令→tτ規(guī)范化時(shí)滯且仍用t表示，則系統(tǒng)（7）變成

dU（t）dt=τ0DΔU（t）+τ0Lμ（Ut）+G（Ut，μ）

（18）

其中

G（Ut，μ）=μDΔU（t）+μLμ（Ut）+（τ0+μ）Fμ（Ut），

Lμ（Ut）=r（u）1-2S*+I*KS-βI*1+αI*Sτ

-r（u）KS*I-βS*（1+αI*）2Iτ

βI*1+αI*Sτ-（δ+γ）I+βS*（1+αI*）2Iτ，

Fμ（Ut）=r（u）（S*-（S2+S*2+SI+S*I*）K）

-β（Sτ+S*）（Iτ+I*）1+α（Iτ+I*）

+β（1+αI*）SτI*+βS*Iτ（1+αI*）2

β（Sτ+S*）（Iτ+I*）1+α（Iτ+I*）

-β（1+αI*）SτI*+βS*Iτ（1+αI*）2

-（δ+γ）I*""""" .

對(duì)系統(tǒng)（18）進(jìn)行線性化得

dU（t）dt=τ0DΔU（t）+τ0Lμ（Ut）

（19）

其線性泛函微分方程為

dz（t）dt=τ0Lμ（zt）

（20）

由Riesz表示定理知，存在一個(gè)有界變差2×2的矩陣函數(shù)η（θ，μ）（θ∈-1，0），滿足

（τ0+μ）Lμ（）=∫0-1dη（θ，μ）（θ），（θ）

∈C（-1，0，R2）

（21）

式（21）中：

dη（θ，μ）=（τ0+μ）Eδ（θ）-（τ0+μ）Fδ（θ+1），

其中

E=m-p

0-z， F=-w-q

wq，

且對(duì)于δ（θ）：［-1，0］→（X，X），有

δ（θ）=0，θ∈［-1，0），

1，θ=0

接下來(lái)定義算子A（0）和A*，對(duì)于（θ）∈C1（［-1，0］，R2），定義A（0）為

A（0）（（θ））=

d（θ）dθ，θ∈［-1，0），

∫0-1dη（θ，0）（θ）， θ=0

對(duì)于Ψ（s）=（Ψ1，Ψ2）∈C1（［0，1］，（R2）*），定義

A*（Ψ（s））=-dΨ（s）ds，s∈［-1，0），

∫0-1Ψ（-ξ）dη（θ，0）， s=0

若Ψ（s）∈C1（［0，1］，（R2）*），（θ）∈C1（［-1，0］，R2），定義如下雙線性形式

（Ψ（s），（θ））0=Ψ-（0）（0）-∫0-1∫θ0Ψ-（ξ-θ）dη（0，θ）（ξ）dξ

（22）

其中，A（0）和A*是正規(guī)伴隨算子.

由于±ω0τ0為A（0）和A*的特征值，算子A（0）對(duì)應(yīng)于iω0τ0的特征向量是q（θ），關(guān)于q的定義如下：

q（θ）=（q1，q2）Teiω0τ0θ（θ∈［-1，0］），

（q1，q2）=（1，m-w-iω0p+q）.

算子A*對(duì)應(yīng)于-ω0τ0的特征向量是q*（s），其中

q*（s）=1（q1*，q2*）Teiω0τ0s（s∈［-1，0］），

（q1*，q2*）=（1，1-m+iω0w），

D=1+q22*-（1-2*）（q1τ0we-iω0τ0+q1q2τ0e-iω0τ0）.令Φ=（q（θ），（θ）），Ψ=（q*（s），*（s））T， 則（Φ，Ψ）0=I， 其中I=10

01.

進(jìn)一步，令f0=（f10，f20）T，其中f10=1

0，

f20=0

1，再定義cf0=c1f10+c2f20，其中c=（c1，c2）T∈C2，在希爾伯特空間X上定義內(nèi)積如下：

〈u，v〉=1lπ∫lπ0u1v-1dx+1lπ∫lπ0u2v-2dx

其中，

u=（u1，u2），v=（v1，v2）∈X=C（［0，lπ］，R2），

并且對(duì)∈C（［-1，0］，X）， 有

〈，f0〉=（〈，f10〉，〈，f20〉）T.

線性化系統(tǒng)（19）的中心子空間為PCNC，且

PCN=Φ（Ψ，〈，f0〉）·f0，∈C，

PCNC={（q（θ）z+（θ））·f0：z∈C}，

對(duì)C進(jìn)行空間分解，故C=PCNCPsC，其中PsC表示穩(wěn)定子空間.

從文獻(xiàn)［22］中可知線性化系統(tǒng)（19）的無(wú)窮小生成元AU滿足AUΨ=Ψ·（θ），其中Ψ∈dom（AU），當(dāng)且僅當(dāng)

Ψ·（θ）∈C，Ψ（0）∈dom（Δ），

Ψ·（θ）（0）=τ0ΔΨ（0）+τ0Lμ（Ψ）.

考慮到分支方向與穩(wěn)定性的判定公式都僅相對(duì)于μ=0，因此在系統(tǒng)（18）中設(shè)μ=0，得到中心流形

W（z，）=W20（θ）z22+W11（θ）z+W02（θ）22+…

（23）

其范圍在PsC中，那么系統(tǒng)（18）在中心流形上的流可以寫成

ut=Φ（z（t），（t））T·f0+W（z（t），（t）），

其中

（t）=iω0τ0z（t）+*（0）〈G（Φ（z（t），（t））T·f0

+W（z，），0），f0〉.

則上式可以寫成（t）=iω0τ0z（t）+g（z，）， 其中

g（z，）=*（0）〈G（Φ（z（t），（t））T·f0+W（z，），0），f0〉

=g20z22+g11z+g0222+g21z22+…

（24）

又因?yàn)镚（，0）=τ0（G1，G1）T=τ0（G1

G1），其中

G1=r（u）（S*-（S2+S*2+SI+S*I*）K）-βS*I*1+αI*-βSτIτ（1+αI*）2

+2αβS*（1+αI*）-3Iτ2+2αβ（1+αI*）-3SτIτ2

-6α2βS*（1+αI*）-4Iτ3，

G2=βSτIτ（1+αI*）2-2αβS*（1+αI*）-3Iτ2+βS*I*1+αI*-（δ+γ）I*

-2αβ（1+αI*）-3SτIτ2+6α2βS*（1+αI*）-4Iτ3

（25）

則根據(jù)式（24）和式（25）可以得到

g20=2τ0D1*-βq1q2e-2iω0τ0（1+αI*）2+2αβS*（1+αI*）-3q22e-2iω0τ0

-r（u）K（q12+q1q2）+2τ0D2*βq1q2e-2iω0τ0（1+αI*）2-2αβS*

（1+αI*）-3q22e-2iω0τ0，

g11=τ0D1*-β（q12+1q2）（1+αI*）2+2αβS*（1+αI*）-32q22-r（u）K（2q11+q12+1q2）

+τ0D2*β（q12+1q2）（1+αI*）2-2αβS*（1+αI*）-32q22，

g02=2τ0D1*-β12e2iω0τ0（1+αI*）2+2αβS*（1+αI*）-322e2iω0τ0-r（u）K（12+12）

+2τ0D2*β12e2iω0τ0（1+αI*）2-2αβS*（1+αI*）-322e2iω0τ0，

g21=-2*1τ0Dr（u）K（1lπ∫lπ02q1W111（0）+1W120（0））dx

+*1τ0D1lπ∫lπ02q1W211（0）+2W120（0）+2q2W111（0）dx

-*1τ0Dβ（1+αI*）21lπ∫lπ0（2q1e-iω0τ0W211（-1）

+1eiω0τ0W220（-1））dx

-*1τ0Dβ（1+αI*）21lπ∫lπ0（2eiω0τ0W120（-1）+2q2e-iω0τ0W111（-1））dx

+4*1τ0DαβS*（1+αI*）-31lπ∫lπ0（2q2e-iω0τ0W211（-1）+2eiω0τ0

W220（-1）+2eiω0τ0W220（-1））dx

+4*1τ0Dαβ（1+αI*）-31lπ∫lπ0（1q22e-iω0τ0+2q1q22e-iω0τ0）dx

-12*1τ0Dα2βS*（1+αI*）-4（2q22e-iω0τ0+22q22e-3iω0τ0）

+*2τ0Dβ（1+αI*）21lπ∫lπ0（2q1e-iω0τ0W211（-1）+1eiω0τ0W220（-1））dx

+*2τ0Dβ（1+αI*）21lπ∫lπ0（2eiω0τ0W120（-1）+2q2e-iω0τ0W111（-1））dx

-4*2τ0DαβS*（1+αI*）-31lπ∫lπ0（2q2e-iω0τ0W211（-1）+2eiω0τ0W220（-1））dx

-4*2τ0Dαβ（1+αI*）-31lπ∫lπ0（1q22e-iω0τ0+2q1q22e-iω0τ0）dx

+12*2τ0Dα2βS*（1+αI*）-4（2q22e-iω0τ0+22q22e-3iω0τ0）.

所以為了計(jì)算g21，需要計(jì)算W20（θ）和W11（θ）.

因?yàn)閃（z（t），（t））滿足

W·=AUW+X0G（Φ（z，）T·f0+W（z，），0）

-Φ（Ψ，〈X0G（Φ（z，）T·f0+W（z，），0），

f0〉）0·f0

=AUW+H20z22+H11z+H0222+…

（26）

根據(jù)W·=W（z，）z+W（z，）·，可得

（2iω0τ0-AU）W20（θ）=H20（θ），

-AUW11（θ）=H11（θ），"""""""""""""""""""" -1≤θlt;0

（-2iω0τ0-AU）W02（θ）=H02（θ）

（27）

當(dāng)-1≤θlt;0時(shí)，

-Φ（Ψ，〈X0G（W（z，）+Φ（z，）T·f0，0），f0〉）·f0

=H20（θ）z22+H11（θ）z+H02（θ）22+….

有

H20（θ）=-［g20q（θ）+02（θ）］·f0，

H11（θ）=-［g11q（θ）+11（θ）］·f0，

代入式（27）得

W20（θ）=ig20ω0τ0q（θ）·f0+i203ω0τ0（θ）·f0+

E1e2iω0τ0θ，

W11（θ）=-ig11ω0τ0q（θ）·f0+i11ω0τ0（θ）·f0+E2.

當(dāng)θ=0時(shí)，利用AU的定義和式（27）可得

H20（0）=-［g20q（0）+02（0）］·f0+τ0E12，

H11（0）=-［g11q（0）+11（0）］·f0+τ0E22，

得到E1=E11·E12，E2=E21·E22，

E11=2iω0-r（u）1-2S*+I*K+βI*e-iω0τ01+αI*r（u）S*K+βS*e-iω0τ0（1+αI*）2

-βI*e-iω0τ01+αI*2iω0+δ+γ-βS*e-iω0τ0（1+αI*）2-1，

E12=-βq1q2e-2iω0τ0（1+αI*）2+2αβS*（1+αI*）-3q22e-2iω0τ0-r（u）K（q12+q1q2）

βq1q2e-2iω0τ0（1+αI*）2-2αβS*（1+αI*）-3q22e-2iω0τ0，

E21=r（u）2S*+I*K-1+βI*1+αI*r（u）S*K+βS*（1+αI*）2

-βI*1+αI*δ+γ-βS*（1+αI*）2-1，

E22=-β（q12+1q2）（1+αI*）2+2αβS*（1+αI*）-32q22-r（u）K（2q11+q12+1q2）

β（q12+1q2）（1+αI*）2-2αβS*（1+αI*）-32q22，

即可確定g21.

基于以上分析，可以看到每個(gè)gij都可由以上參數(shù)來(lái)確定，則以下用于判斷Hopf分支方向以及分支周期解性質(zhì)的參數(shù)μ2，β2，T2的值可被計(jì)算出

C1（0）=i2ω0τ00（g11g20-2g112-g0223）+g212，

μ2=-Re（C1（0））Re（λ′（τ00）），β2=2Re（C1（0）），

T2=-Im（C1（0））+μ2Im（λ′（τ00））ω0τ00.

由此得到定理2.

定理2 對(duì)系統(tǒng)（5），

（1）μ2決定Hopf分支方向：若μ2gt;0（μ2lt;0），則分支周期解產(chǎn)生于τgt;τ00（τlt;τ00）的一側(cè).

（2）β2決定Hopf分支周期解的穩(wěn)定性：若β2lt;0，則分支周期解穩(wěn)定; 若β2gt;0，則分支周期解不穩(wěn)定.

（3）T2決定Hopf分支周期解周期的變化：若T2gt;0，則周期變大; 若T2lt;0，則周期變小.

3 數(shù)值模擬

在本節(jié)中，將選取兩組滿足系統(tǒng)正平衡點(diǎn)存在且能使系統(tǒng)（6）產(chǎn)生一對(duì)純虛根的數(shù)據(jù)進(jìn)行數(shù)值模擬來(lái)驗(yàn)證本文的理論分析.

對(duì)于系統(tǒng)（6），r0=0.95，K=2，β0=1.4，α=0.01，δ+γ=0.95，dS=1，dI=0.1.經(jīng)計(jì)算可知ω0=0.677 2， 地方病平衡點(diǎn)E*（S*，I*）=（0.714，0.338）且滿足（δ+γ）lt;β（u，v）lt;r（u）+（δ+γ）.此時(shí)有τ00=0.4915.由定理1和2可知，

當(dāng)τlt;τ00時(shí)， 系統(tǒng)在地方病平衡點(diǎn)E*（S*，I*）處是趨于局部漸近穩(wěn)定的， 如圖1所示;當(dāng)τ經(jīng)過τ00=0.491 5時(shí)，地方病平衡點(diǎn)E*（S*，I*）會(huì)失去它的穩(wěn)定性并產(chǎn)生Hopf分支，如圖2所示.由定理1和第1節(jié)公式推導(dǎo)得β2=-6.4×10-3lt;0，μ2=1.46×10-2gt;0，T=2×10-1gt;0，由此可知當(dāng)τ經(jīng)過τ00=0.491 5時(shí)，產(chǎn)生的Hopf分支為超臨界，并且分支周期解是穩(wěn)定的，周期增加.

取τ=0.3lt;τ00和初值S（x，t）=0.714+0.5 cos5x，I（x，t）=0.338+0.5cos5x，由定理1可知系統(tǒng)在E*（S*，I*）處是局部漸近穩(wěn)定，如圖1所示.

取τ=0.6gt;τ00，由定理1和2可知系統(tǒng)在地方病平衡點(diǎn)E*（S*，I*）處產(chǎn)生Hopf分支.當(dāng)初值取S（x，t）=0.714+0.5cos5x，I（x，t）=0.338+0.5cos5x時(shí)，系統(tǒng)產(chǎn)生空間齊次分支周期解，如圖2所示.

4 結(jié)論

本文提出了一個(gè)具有進(jìn)化效應(yīng)的時(shí)滯反應(yīng)擴(kuò)散傳染病模型，研究了時(shí)滯對(duì)該模型動(dòng)力學(xué)行為的影響.通過分析發(fā)現(xiàn)，當(dāng)時(shí)滯參數(shù)經(jīng)過某一臨界值時(shí)，會(huì)改變系統(tǒng)的穩(wěn)定性，此時(shí)系統(tǒng)在正平衡點(diǎn)附近出現(xiàn)周期振蕩模式并產(chǎn)生Hopf分支.同時(shí)通過Matlab給出具體的數(shù)值實(shí)例，可以看出當(dāng)時(shí)滯參數(shù)小于某個(gè)臨界值時(shí)，易感人群和感染人群的密度會(huì)在特定區(qū)域內(nèi)趨于一個(gè)正平衡態(tài)，保持穩(wěn)定; 而當(dāng)時(shí)滯參數(shù)超過該臨界值時(shí)，兩人群的密度變化會(huì)呈現(xiàn)周期振蕩模式.

因此，時(shí)滯效應(yīng)是影響傳染病系統(tǒng)人群密度變化的一個(gè)重要因素.

參考文獻(xiàn)

［1］ Kermack W O，Mckendrick A G.A contribution to the mathematical theory of epidemics［J］.Bulletin of Mathematical Biology，1991，53（1-2）：57-87.

［2］ Kermack W O，Mckendrick A G.Contributions to the mathematical theory of epidemics II：The problem of endemicity［J］.Proceedings of the Royal Society of London（Series A），1932，138（834）：55-83.

［3］ Ruan S G，Wang W D.Dynamical behavior of an epidemic model with a nonlinear incidence rate［J］.Journal of Differential Equations，2003，188（1）：135-163.

［4］ Wang J J，Zhang J Z，Jin Z.Analysis of an SIR model with bilinear incidence rate［J］.Nonlinear Analysis Real World""" Applications，2010，11（4）：2 390-2 402.

［5］ Xiao Y N，Chen L S，Bosch F V，et al.Dynamical behavior for a stage-structured SIR infectious diseasemodel［J］.Nonlinear Analysis：Real World Applications，2002，3（2）：175-190.

［6］ Huang G，Takeuchi Y.Global analysis on delay epidemiological dynamic models with nonlinear incidence［J］.Journal of Mathematical Biology，2011，63（1）：125-139.

［7］ Jeschke J M.Density-dependent effects of prey defenses and predator offenses［J］.Journal of Theoretical Biology，2006，242（4）：900-907.

［8］ Zhao Q Y，Liu S T，Niu X L.Dynamic behavior analysis of a diffusive plankton model with defensive and offensive effects［J］.Chaos，Solitions amp;Fractals，2019，129：94-102.

［9］ Velzen E V，Gaedke U.Reversed predator-prey cycles are driven by the amplitude of prey oscillations［J］.Ecology and Evolution，2018，8（12）：6 317-6 329.

［10］ Velzen E V，Gaedke U.Disentangling eco-evolutionary dynamics of predator-prey coevolution：The case of antiphase cycles［J］.Scientific Reports，2017，7（1）：11.

［11］ 楊 洪，馬秋敏，盛江林，等.具有進(jìn)化效應(yīng)的SIR模型的穩(wěn)定性和Hopf分支分析［J］.高校應(yīng)用數(shù)學(xué)學(xué)報(bào)：A輯，2023，38（1）：27-36.

［12］ Anderson R M，May R M.Population biology of infectiousdiseases［J］.Nature，1979，180：361-367.

［13］ Abdelilah K，Abdelhadi A，Hamad T A.A comparison of delayed SIR and SEIR epidemic models［J］.Nonlinear Analysis Modelling and Control，2011，16（2）：181-190.

［14］ Yan P，Liu S Q.SEIR epidemic model with delay［J］.Anziam Journal，2006，48（1）：119-134.

［15］ Kuang Y.Delay differential equations：With applications in population dynamics［M］.New York：Academic Press，1993.

［16］ Neuhauser C.Mathematical challenges in spatial ecology［J］.Notices of the American Mathematical Society，2001，48（11）：1 304-1 314.

［17］ Murray J D.Mathematicalbiology II：Spatial models and biomedical applications［M］.Berlin：Springer，2003.

［18］ Rass L，Radcliffe J.Spatial deterministic epidemics［M］.Rhode Island：American Mathematical Society，2003.

［19］ Ruan S G，Wu J H.Modeling spatial spread of communicable diseases involving animal hosts［M］.New Jersey：Spatial Ecology，2009.

［20］ Wu J H.Theory and applications of partial functional differential equation［M］.Berlin：Springer，1996.

［21］ Lin X.Centre manifolds for partial differential equations with delays［J］.Proceedings of the Royal Society of Edinburgh：Section A Mathematics，1992，122（3-4）：237-254.

［22］ Jiang W H，Wang H B，Cao X.Turing instability and turing-hopf bifurcation in diffusive schnakenberg systems with gene expression time delay［J］.Journal of Dynamics and Differential Equations，2018，31（4）：2 223-2 247.

【責(zé)任編輯：陳 佳】