張琬婧，林支桂

(揚州大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，江蘇 揚州 225002)

近期H1N1[1]、H7N9[2]、登革熱[3]等傳染病反復(fù)發(fā)作，新型冠狀病毒肺炎[4]不斷流行，傳染病的傳播引起流行病學(xué)專家以及社會廣泛關(guān)注。傳染病受季節(jié)、溫度等因素影響，其傳播范圍隨時間變化，意味著受傳染病影響的區(qū)域與時間有關(guān)?？紤]到種群在空間上的遷徙，在模型中研究種群密度的時空演化能更好地描述傳染病的變化規(guī)律，故而在常微分方程的基礎(chǔ)上引入擴散作用，從而得到反應(yīng)擴散方程。此類方程描繪了生物種群隨時空變化而產(chǎn)生的一系列規(guī)律[5]。但是引入擴散作用后，原本模型的一系列穩(wěn)定性質(zhì)就會被打亂，其中Turing不穩(wěn)定就是一個例子：由均一的平衡狀態(tài)變?yōu)椴痪粻顟B(tài)，這一工作得到了人們的極大關(guān)注[6-9]。由此，為理解在增長的棲息地區(qū)域中，擴散對物種生存的影響，本文在增長區(qū)域下討論一類寄生蟲-宿主傳染病模型的Turing不穩(wěn)定。

1 擴散的寄生蟲-宿主傳染病模型

2004年，在Kermack等[10]建立的經(jīng)典倉室模型基礎(chǔ)上，Berezovsky等[11]提出1個包含可變?nèi)丝?、接觸傳播、因病死亡等因素的寄生蟲-宿主模型。

(1)

式中：N表示物種總數(shù)，它包含易感人群S和染病人群I，且滿足N=S+I；r表示易感人群S的內(nèi)稟增長率；K表示易感人群S的承載能力；β表示染病人群I的接觸傳播率；μ表示自然死亡率；d表示病死率；m表示非染病者的人均移民率。

Wang等[12]將模型(1)進行無量綱化，并引入擴散

(2)

其對應(yīng)的常微分方程的正平衡解所對應(yīng)的Jacobi矩陣為

下面討論系統(tǒng)(2)解的Turing失穩(wěn)條件。

定理1當1

(d1a22+d2a11)2-4d1d2(a11a22-a12a21)>0，

且存在某個δi>0，使得不等式0

{φij,j=1,2,…,dimE(δi)}是E(δi)的一組標準正交基，

μ是Xi上的特征值當且僅當μ是矩陣

的特征值。此時特征方程為

顯然，要證明E*對于系統(tǒng)(2)不穩(wěn)定，只需要det(Ai)<0，也就是特征方程至少有1個正解。因為δi>0，所以不穩(wěn)定第一個必要條件就是d1a22+d2a11>0，這等價于a11>-d1a22/d2>0。

(d1a22+d2a11)2-4d1d2(a11a22-a12a21)>0。

綜上所述，要使系統(tǒng)(2)在平衡解E*處Turing不穩(wěn)定，那么Rd、R0、ν必須滿足

(d1a22+d2a11)2-4d1d2(a11a22-a12a21)>0。

事實上，此情形還不能保證一定存在δi>0，使得

(3)

所以直接計算

式中M=detJ=a11a22-a12a21。

在滿足上述條件的情況下，方程h(δi)=0存在2個正解：

從而有結(jié)論：如果存在δi>0，使不等式0

定理2假設(shè)E*是系統(tǒng)(2)的染病平衡點，如果Ω=(0,lπ)，滿足

那么系統(tǒng)(2)在平衡解E*處Turing不穩(wěn)定。

2 增長區(qū)域上的寄生蟲-宿主傳染病模型

第1章通過線性化方法給出了固定區(qū)域下寄生蟲-宿主傳染病模型的Turing不穩(wěn)定的條件，結(jié)果表明，擴散會導(dǎo)致此模型的Turing不穩(wěn)定。本章將考慮增長區(qū)域上具Neumann邊界條件的寄生蟲-宿主模型。記

運用Plaza等[13]建立的框架，將區(qū)域的增長率和曲率融入傳染病動力學(xué)模型，即設(shè)二維平面X嵌入三維空間R3中，表示為

X(ξ,η,t)=(x(ξ,η,t),y(ξ,η,t),z(ξ,η,t))，

并定義

(4)

(5)

式中：S0(X(0))和I0(X(0))是正有界函數(shù)；Ω(0)為初始區(qū)域。由此得到增長區(qū)域下具Neumann邊界條件下的寄生蟲-宿主傳染病模型。

為刻畫區(qū)域演化率的影響，考慮棲息地的演化是各向同性的。首先探索限制在平面區(qū)域上增長時的情形，數(shù)學(xué)上表達式為

X(ξ,η,t)=ρ(t)Y(ξ,η)=ρ(t)[ξ,η,0]T，

式中ρ(t)是增長函數(shù)，滿足ρ(t)在[0，+∞)上連續(xù)可微，

經(jīng)計算，得

由此，式(5)轉(zhuǎn)化為平面增長區(qū)域下的寄生蟲-宿主傳染病擴散模型

(6)

3 增長區(qū)域上的Turing不穩(wěn)定

(7)

從而式(6)變?yōu)?/p>

(8)

式(8)對應(yīng)的動力系統(tǒng)為

(9)

證明為了得到具有每種分支類型(Hopf、Turing和Turing-Hopf)下的參數(shù)值的條件，根據(jù)文獻[14]，先做下列計算。式(7)對應(yīng)的ODE系統(tǒng)的零斜率線為

(10)

I1和I2相交所滿足的等式為

(11)

因此，給出條件

(12)

從而可以保證式(11)中SΔ>0，IΔ>0。故EΔ=(SΔ,IΔ)是正平衡解。

其次，記(SΔ,IΔ)處的Jacobi矩陣為J(EΔ)，這里

(13)

式中：

考慮其特征方程

(14)

經(jīng)計算，

類似定理1的證明，通過線性化得特征矩陣，討論特征方程，給出其至少存在1個正根的必要條件，即得系統(tǒng)產(chǎn)生Turing不穩(wěn)定的必要條件。

定理4式(8)產(chǎn)生Turing不穩(wěn)定的必要條件為

fS+gI-4k<0，

fSgI-fIgS-2k(fS+gI)+4k2>0，

(d2fS+d1gI)-2k(d1+d2)>0，

把f、g的表達式代入計算可得

1+2k

(15)

(16)

(17)

(18)

由此，通過分析發(fā)生Turing不穩(wěn)定的必要條件，可以發(fā)現(xiàn)：擴散系數(shù)d2越大，模型更容易出現(xiàn)斑圖模式；而區(qū)域增長快，即k大時，上述必要條件不容易滿足，所以對Turing斑圖起破壞作用。

4 數(shù)值模擬及解釋

本章將通過數(shù)值模擬來研究固定區(qū)域上Turing不穩(wěn)定在此類寄生蟲-宿主模型中的動力學(xué)行為。在實際流行病學(xué)中，空間上的擴散影響物種初始活動空間范圍[15]。為研究模型非平凡解的動態(tài)分布，選取參數(shù)R0=1.3，Rd=1.7，ν=0.3，d1=0.1。

對時間和空間都采用有限差分方法，且取時間步長Δt=0.01，空間步長Δx=1.5，最終時間t=1 000，初始條件是種群初值的隨機分布。

由圖1可以得出，當擴散系數(shù)d2越小，即(d1a22+d2a11)2-4d1d2(a11a22-a12a21)<0時，種群最后會收斂到一個平衡解，不發(fā)生Turing不穩(wěn)定；當擴散系數(shù)d2較大時，Turing不穩(wěn)定的充分條件滿足，因此種群在不同位置的值不同，于是常數(shù)解不穩(wěn)定。圖2則是產(chǎn)生Turing不穩(wěn)定時種群S和I所產(chǎn)生的斑圖。

圖1 不同擴散系數(shù)下種群S在位置(25，25)、(50，50)、(75，75)時的密度Fig. 1 Densities of population S at locations (25, 25), (50, 50) and (75, 75) with different diffusion coefficients

圖2 d2=5時種群S和I所形成的斑圖Fig. 2 Spatial patterns of population S and I when d2=5