周期環(huán)境下具有潛伏期和階段結(jié)構(gòu)的時滯布魯氏菌病模型研究

吳 浩,桑 瑞,張 龍,李洪利,滕志東

(新疆大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,新疆烏魯木齊 830046)

§1 引言

布魯氏菌病(brucellosis)是人類和動物共患的一種具有嚴(yán)重危害的傳染病,它是全世界范圍內(nèi)常見的由細(xì)菌引起的疾病[1-4].近年來,許多學(xué)者通過構(gòu)造數(shù)學(xué)模型來研究布魯氏菌病的傳播規(guī)律,并取得了一些重要的研究成果.Hou[5]等人在2013年提出了一類自治常微分方程模型用來描述布魯氏菌病的傳播動態(tài),作者考慮了布魯氏菌病直接傳播和間接傳播兩種疾病傳播方式,從而建立了SVEIB模型,確定了基本再生數(shù)R0,分析了無病平衡點(diǎn)和地方病平衡點(diǎn)的全局穩(wěn)定性.

2014年Sun[6]等人考慮了一類具有遷移的自治SEIVW常微分模型

其中他們將動物種群分為四類:易感者類,暴露者類,感染者類,接種者類,以及一類環(huán)境中的布魯氏菌類,分別用S(t),E(t),I(t),V(t)和W(t),表示它們在t時刻的種群密度.模型中其他參數(shù)的意義分別為:A,遷入率;d和α,自然死亡率和因病死亡率;β,β1,k,動物之間的直接接觸的疾病傳播率,動物與布魯氏菌接觸的疾病傳播率和病原體的脫落率;r和δ,疫苗接種率和免疫下降率;nτ和?,環(huán)境的去污率和病原體的衰減率;b(S(t)+V(t)),出生率.他們確定了該模型的基本再生數(shù)R0,并建立完全由R0確定的全局動力學(xué),即如果R0<1,模型(1)的無病平衡點(diǎn)全局漸近穩(wěn)定;而R0>1,模型(1)存在一個全局漸近穩(wěn)定的地方病平衡點(diǎn).2016年,Hou[7]等人提出了具有階段結(jié)構(gòu)的布魯氏菌病模型,性未成熟的動物幾乎不會通過直接感染或間接感染的方式感染布魯氏菌病[8-9],故他們考慮將易感者按年齡分類即分為幼年動物和成年動物.他們確定了模型的基本再生數(shù)R0,通過對基本再生數(shù)關(guān)于一些參數(shù)進(jìn)行敏感性分析,發(fā)現(xiàn)了羊的出生率,接種率,以及具有感染性的羊的清除率在布魯氏菌病的傳播中起著重要作用.

潛伏期是布魯氏菌病研究中考慮的重要因素之一.許多受感染的動物都會表現(xiàn)出較長的潛伏期,在此期間,動物的血清學(xué)表現(xiàn)為非陽性,即潛伏期的感染動物幾乎不可能通過檢測被發(fā)現(xiàn)[10-12].考慮布魯氏菌病傳播的感染潛伏期,Hou[11]等人既引入感染潛伏期時滯建立SVIP模型,也通過引入暴露類E建立SVEIP模型,利用李雅普諾夫函數(shù)和泛函分別研究了兩個模型平衡點(diǎn)的穩(wěn)定性，發(fā)現(xiàn)兩個模型的動力學(xué)性質(zhì)完全依賴于各自的基本再生數(shù).2019 年Hou[12]等人建立了一類布魯氏菌傳染病模型,該模型包含了易感個體的不同階段和感染時滯,結(jié)果表明,平衡點(diǎn)的全局動力學(xué)依賴于基本再生數(shù)R0:若R0≤1，動物布魯氏菌病將最終滅絕;若R0>1,動物布魯氏菌病是持久的,最終趨向于穩(wěn)定的地方病平衡點(diǎn).

值得注意的是上面的工作都是自治環(huán)境下研究的,而布魯氏菌病傳播過程中的季節(jié)性,特別是氣候因素,如溫度,濕度和紫外線照射強(qiáng)度會顯著影響布魯氏菌在環(huán)境中的生存[13-15].因此,研究布魯氏菌病與環(huán)境季節(jié)性的關(guān)系,是其動力學(xué)研究的基礎(chǔ),而這方面的研究較少[16-18],2017年P(guān).O.Lolika[16]等人研究了一類季節(jié)性的SHEIB布魯氏菌病模型即模型(2),

他們確定了與時間周期模型相關(guān)的基本再生數(shù),以其作為閾值,討論了模型(2)的全局動態(tài),發(fā)現(xiàn)季節(jié)性在塑造布魯氏菌病的長期動態(tài)中起著重要的作用,從而影響其最優(yōu)控制策略的設(shè)計(jì).Lou[17]等人建立了一類具有周期性傳播率的SEIV模型用于研究布魯氏菌病的傳播,研究表明,提高感染羊或牛的屠宰率，提高易感羊或牛的接種率，降低疫苗的損失率是控制布魯氏菌病流行的有效策略.2021年Yang[18]等人提出了一類具有非局部傳播和空間擴(kuò)散的季節(jié)性布魯氏菌病流行模型,定義了基本再生數(shù)R0,結(jié)果表明,擴(kuò)大有效感染半徑增加了布魯氏菌病傳播的風(fēng)險.建議牧民應(yīng)改善動物環(huán)境衛(wèi)生,及時隔離感染動物,從而達(dá)到控制布魯氏菌病的流行.

基于前人的工作,本文基于以下幾個方面建模,將易感者按年齡分類即分為幼年動物和成年動物;引入布魯氏菌病傳播的感染潛伏期;由于布魯氏菌病可導(dǎo)致流產(chǎn)且大大降低了幼畜的存活率,所以假設(shè)感染的動物不能孕育后代,動物接種和感染動物的捕殺作為控制布魯氏菌病的措施.本文的目的是在周期模型(2)的基礎(chǔ)上,將上述因素引入周期模型,使其更符合布魯氏菌病傳播的實(shí)際情況，從而提出一類周期環(huán)境下的一般布魯氏菌病模型.然后,對建立的周期模型進(jìn)行動力學(xué)分析,重點(diǎn)關(guān)注其與基本再生數(shù)相關(guān)的閾值動態(tài)特征.

§2 模型的建立和介紹

為了更符合現(xiàn)實(shí)的布魯氏菌病傳播機(jī)制和進(jìn)一步研究其動力學(xué)行為,本文建立了如下一類周期環(huán)境下具有潛伏期及階段結(jié)構(gòu)的時滯布魯氏菌病模型

將動物種群分為四類:幼年易感者類,成年易感者類,接種者類,感染者類,以及一類環(huán)境中的布魯氏菌類,其在t時刻的數(shù)量和濃度分別用S1(t),S2(t),H(t),I(t)和B(t)表示.遷入率用A表示,幼年動物的自然死亡率,成年動物的自然死亡率和成年動物的因病死亡率分別用μ1,μ2,c表示.幼年動物向成熟個體的轉(zhuǎn)化率用σ1表示.動物之間的直接接觸的疾病傳播率,動物與布魯氏菌接觸的疾病傳播率和病原體的脫落率分別用β1(t),β2(t),β3(t)表示.疫苗接種率和免疫下降率分別用l和k表示,修正因子,感染潛伏期和感染動物的捕殺率分別用γ,τ,p表示.環(huán)境的去污率和病原體的衰減率分別用δ和d(t)表示.出生率用表示.

§3 模型的基本性質(zhì)

引理3.1模型(3)滿足初始條件(4)有唯一的,非負(fù)解

在[0,∞)上有定義,其中φ ∈C+.

證對于任意的φ ∈C+,定義

由于f(t,φ)在(t,φ)∈R×C+上連續(xù),且f(t,φ)關(guān)于φ在R×C+所有緊集上是Lipschitz的,所以由文獻(xiàn)[19]定理2.2.1,定理2.2.3即泛函微分方程解的存在唯一性可知:模型(3)滿足u0=φ有唯一解u(t,φ)存在于它的最大存在區(qū)間[0,σφ)上.

若φi(0)=0,i ∈{1,2,3,4,5},則fi(t,φ)≥0,故由文獻(xiàn)[20]定理5.2.1可知:模型(3)滿足u0=φ的唯一解u(t,φ)有:ut(φ)∈C+成立,其中t ∈[0,σφ).

令M(t)=S1(t)+S2(t)+H(t)+eμτI(t+τ),則

其中m=min{μ1,μ2}.進(jìn)一步可得:M(t)在[0,σφ)上有界,因此S1(t),S2(t),H(t),I(t)在[0,σφ)上有界,從而B(t)有界.根據(jù)文獻(xiàn)[19]定理2.3.1即泛函微分方程延拓定理可知:σφ=∞,即S1(t),S2(t),H(t),B(t),I(t)在[0,∞)上有定義且非負(fù).

對于模型(3),定義了模型的解半流P(t) :C+→C+,即P(t)φ=ut(φ),其中φ ∈C+.與文獻(xiàn)[21]引理3.5相似的證明可得:

引理3.2模型(3)產(chǎn)生一個周期半流(periodic semiflow)P=P(ω),且有一個全局吸性子(global attractor).

證顯然P0=I,Pt(φ)連續(xù).故只需證明Pt+ω=Pt °Pω.定義x(t)=u(t+ω,φ),下面證明x(t)=u(t,uω(φ)).

所以x1是模型(3)的解,同理用x(t)=u(t+ω,φ)和βj(j=1,2,3)的周期性可得:x2,x3,x4,x5是模型(3)的解,故x(t)是模型(3)的一個解.

顯然xi(θ)=ui(θ+ω,φ)成立,其中θ ∈[?τ,0],i=1,2,3,4,5.

另一方面,令y(t)=u(t,uω(φ)),所以y(t)也是模型(3)的一個解.而yi(θ)=ui(θ,uω(φ))=ui(θ+ω,φ).所以x(t)和y(t)都是模型(3)滿足相同初始條件(4)的解.由解的唯一性可知:x(t)=y(t),即u(t+ω,φ)=u(t,uω(φ))成立,其中t ≥0.

其中t ≥0,θ ∈[?τ,0].同理有ut°uω(φi)(θ)=ut+ω(φi)(θ)成立,其中θ ∈[?τ,0],i=2,3,4,5.所以ut°uω(φ)=ut+ω(φ).故Pt+ω(φ)=Pt °Pw(φ)成立,其中φ ∈C+,t ≥0成立.綜上所述,{P(t)}t≥0是C+上的ω周期半流.

由前面的引理3.1的證明易得

所以模型(3)的解S1(t),S2(t),H(t),I(t)都是最終有界的,從而B(t)最終有界,即解半流P(t) :C+→C+是點(diǎn)耗散的(point dissipative).由文獻(xiàn)[19]定理3.6.1可知,對于任意t>τ,P(t)是緊的,進(jìn)而通過文獻(xiàn)[22]定理3.1.1可知P=P(ω)在C+有全局吸引子.

§4 基本再生數(shù)

接下來根據(jù)Zhao等人在文獻(xiàn)[23]發(fā)展的理論來推導(dǎo)模型(3)的基本再生數(shù).只需要考慮感染相關(guān)的倉室,即I倉室和B倉室,在E0處對其線性化得到

令Z=C([?τ,0],R2),Z+=C([?τ,0],R2+),(Z,Z+)是一個有序的Banach空間,有最大模范數(shù)和正錐Z+.對于函數(shù)z(·)=(z1(·),z2(·))∈C([?τ,∞),R2),定義泛函zt ∈Z,zt(θ)=(z1(t+θ),z2(t+θ)),其中θ ∈[?τ,0].故

在t時刻的新感染個體線性依賴于時間間隔[t ?τ,t]的感染個體,可以通過F(t)zt來描述.個體在受感染倉室的內(nèi)部進(jìn)化(例如:自然死亡,因病死亡等)由線性常微分模型

描述.令Φ(t,s),t ≥s是其相關(guān)聯(lián)的進(jìn)化矩陣即對于任意給定s ∈R,2×2矩陣Φ(t,s)滿足

其中I是2×2的單位陣.

易得

考慮到Φ(t,s)指數(shù)增長的界:

w(Φ):=inf{w:存在M ≥1,使得||Φ(t+s,s)||≤Mewt,其中s ∈R,t ≥0}.令d(t)+δ,則

(H1):所有算子F(t):Z →R2是正的(F(t)Z+?);

(H2):周期矩陣?V(t)是合作的(cooperative)且w(Φ)<0.

令Cω是所有從R映到R2的ω周期函數(shù)的Banach空間,具有最大模范數(shù)和正錐={u ∈Cω:u(t)≥0,其中t ∈R}.根據(jù)文獻(xiàn)[22]的理論,定義Cω上的線性算子L,

對于任意給定的s ≥0,vt?s是在t?s時刻感染個體的初始分布,F(t?s)vt?s是在t?s時刻新感染個體的分布,它們是由在區(qū)間[t ?s ?τ,t ?s]的感染個體產(chǎn)生的.Φ(t,t ?s)F(t ?s)v(t ?s+·)是在t ?s時刻新感染且在t時刻仍保留在受感染倉室的感染個體分布.

即t時刻新感染的累積分布,這些新感染個體由t時刻之前所有感染個體產(chǎn)生.故定義基本再生數(shù)為

即基本再生數(shù)是算子L的譜半徑(spectral radius).

令P1(t)是I(t),B(t)在無病平衡點(diǎn)E0處的線性化方程的解映射,即P1(t)φ=zt(φ),t ≥0,其中z(t,φ)是模型(6)滿足z0=φ ∈Z的唯一解,P1=P1(ω)是模型(6)相關(guān)聯(lián)的龐加萊映射(Poincaré map),令r(P1)是P1的譜半徑.根據(jù)文獻(xiàn)[23]定理2.1有

引理4.1R0?1和r(P1)?1有相同的符號.

考慮下面的相空間

給定一個連續(xù)函數(shù)α:[?τ,∞]×R→R2,定義αt ∈W,

與前面的方法一樣,可以推出模型(6)在[0,∞)有唯一滿足α0=φ非負(fù)的,最終有界的解α(t,φ),其中φ ∈W+.令P2:=P2(ω)是模型(6)的相關(guān)聯(lián)的龐加萊映射,r(P2)是P2的譜半徑.

與文獻(xiàn)[24]引理3.8相似的證明可得

引理4.2r(P1)=r(P2).

§5 布魯氏菌病的滅絕與持久

首先給出了一些準(zhǔn)備工作,下面證明周期半流P2(t)是最終強(qiáng)正的(eventually strongly positive).

引理5.1若α(t,φ)是模型(6)滿足α0=φ的解,則αi(t,φ)>0成立,其中t >2τ,i=1,2且φ ∈W+{0},進(jìn)一步有P2(t)φ ?0成立,其中t>3τ.

證令α(t,φ)=(α1(t,φ),α2(t,φ))是模型(6)滿足α0=φ的解,其中φ=(φ1,φ2)∈W+{0}.為了簡化符號,定義(α1(t),α2(t))=(α1(t,φ),α2(t,φ)).不失一般性,假設(shè)φ1>0,下面證明:存在t0∈[0,τ],使得α2(t0)>0.反證.假設(shè)α2(t)=0成立,其中t ∈[0,τ].

由模型(6)的第2個方程,

綜上所述,αt(φ)?0成立,其中t>3τ.故P2(t)φ ?0成立,其中t>3τ.

接下來,給出本節(jié)的兩個重要結(jié)論:布魯氏菌病的滅絕和一致持久(uniformly persistent).

定理5.3當(dāng)R0<1時,模型(3)的無病平衡點(diǎn)E0在C+上全局吸引.

證考慮R0<1,令P?是下面輔助模型的龐加萊映射.

所以存在一個充分大的整數(shù)n1>0且滿足n1ω ≥τ,使得

對于任意的t ≥n1ω,有

令X0:={φ=(φ1,φ2,φ3,φ4,φ5)∈C+:φ4(0)>0 且φ5(0)>0};

?X0:={φ=(φ1,φ2,φ3,φ4,φ5)∈C+:φ4(0)=0 或φ5(0)=0}.

定理5.4當(dāng)R0>1時,則模型(3)有一個正的ω周期解,且模型(3)是一致持久的,即存在一個實(shí)數(shù)η >0,使得對于模型(3)任意正解

(S1(t),S2(t),H(t),I(t),B(t))=(S1(t,φ),S2(t,φ),H(t,φ),I(t,φ),B(t,φ)),

其中初值滿足φ ∈X0,有

證考慮R0>1,則r(P2)>1.Q(t)是模型(3)在C+上的解映射,即Q(t)φ=ut(φ),t ≥0,Q:=Q(ω)是模型(3)相關(guān)聯(lián)的龐加萊映射,顯然,Q(t)X0?X0成立,其中t ≥0.下面證明Q關(guān)于(X0,?X0)一致持久.

令Pδ是下面擾動輔助模型的龐加萊映射.

故對于任意的t ≥n2ω+τ,由模型(3)的第4,第5個方程,即I(t,φ),B(t,φ)滿足

由于Q(t)X0?X0,其中t ≥0,所以可以選擇一個充分小的k >0,使得

故根據(jù)泛函微分方程的比較原理可知

定義M?:={φ ∈?X0:Qn(φ)∈?X0,其中n ≥0}.Qn(φ)∈?X0成立,其中φ ∈M?且n ≥0.所以對于任意的n ≥0,I(nω,φ)=0,B(nω,φ)=0至少有一個成立.

由前面引理3.1的證明可知:若存在t0≥0,使得I(t0)>0,則I(t)>0成立,其中t ≥t0;若存在t0≥0,使得B(t0)>0,則B(t)>0成立,其中t ≥t0.故I(t,φ)≡0成立,或者B(t,φ)≡0成立,其中t ≥0.

若I(t,φ)≡0,其中t ≥0.由模型(3)的第5個方程,顯然有B(t)→0,t →∞,進(jìn)而有S1(t)→0,S2(t)→0,H(t)→0,t →∞;若存在t0≥0使得I(t0,φ)>0,進(jìn)而有

所以B(t,φ)≡0,其中t ≥t0,故S1(t)→0,S2(t)→0,H(t)→0,t →∞.

所以ω(φ)=M1,其中φ ∈M?.即M1關(guān)于Q在?X0上不能形成環(huán).由一致持久性的無環(huán)性定理即文獻(xiàn)[22]定理1.3.1可知:Q:C+→C+關(guān)于(X0,?X0)一致持久.又由于Qn是緊的,其中整數(shù)n滿足nω >τ,根據(jù)文獻(xiàn)[22]定理1.3.6可知

§6 數(shù)值模擬

下面通過數(shù)值例子來模擬布魯氏菌病的傳播過程,取τ=0.2;A=1000;b=0.8;?=0.4;μ1=0.15;σ1=0.75;μ2=0.25;k=0.75;γ=0.45;c=0.102;δ=0.0275.首先,圖1顯示了布魯氏菌病一致持久.圖2顯示了布魯氏菌病滅絕,從圖2中可以看出，從不同初始值開始的解最終趨于無病平衡點(diǎn).

圖1 τ=0.2時,布魯氏菌病一致持久,(a)各倉室的時間序列圖;(b) S2(t),I(t),B(t)倉室的三維相圖;(c)S2(t)倉室的時間序列圖;(d) I(t) 倉室的時間序列圖;l=0.29;p=0. 200025;

圖2 τ=0.2時,布魯氏菌病滅絕時,各倉室的時間序列圖,l=0.89;p=0. 700025;

接著,研究時滯τ以及周期性對模型(3)的影響.從圖3(a)中可以看到,隨著τ的增加,I(t)的最大值增加,最小值減小,疾病的波動變大.從圖3(b)中可以看到,無周期時I(t)最后趨于一條直線即地方病平衡點(diǎn),在模擬的情況中ω=12的峰值最大,波動最大.它們波動的平均值也是有差異的,也就是說忽略時滯和周期性,疾病的風(fēng)險可能被低估或高估.對于短期的研究自治模型適用,而對于布魯氏菌病長期的研究周期模型更為適用,本文考慮布魯氏菌病的潛伏期和季節(jié)性可能與現(xiàn)實(shí)更吻合.

圖3 (a) 其他參數(shù)不變,分別取τ=0.01,τ=50,τ=200, I(t)的時間序列圖;(b) 其他參數(shù)不變,分別取ω=365, ω=150,ω=12以及無周期, I(t) 的時間序列圖

本文主要闡述接種和感染動物的捕殺率對于布魯氏菌病的影響,這兩種措施是容易實(shí)現(xiàn)且有效的.從圖4中可以看到隨著接種和感染動物的捕殺的增加,布魯氏菌病可以得到有效的控制.通過控制接種和感染動物的捕殺率可以使得原來一致持久的布魯氏菌病滅絕,當(dāng)然現(xiàn)實(shí)情況這兩個參數(shù)值增大到一定程度可能難以甚至無法實(shí)現(xiàn),在本文研究的基礎(chǔ)上還可以考慮對其進(jìn)行最優(yōu)控制,在本文暫未考慮,可以放在以后的工作中繼續(xù)討論.

圖4 (a) 其他參數(shù)不變,分別取l=0.19,l=0.29,l=0.59, I(t)的時間序列圖;(b) 其他參數(shù)不變,分別取p=0.000025,p=0.100025,p=0.200025以及p=0.200025, I(t) 的時間序列圖

其次,如果不考慮動物的階段結(jié)構(gòu),有以下模型:

在相同的參數(shù)條件下,分別模擬有無階段結(jié)構(gòu)模型.從圖5可以看出,無階段結(jié)構(gòu)時疾病滅絕,考慮階段結(jié)構(gòu)時疾病一致持久,因此如果忽略動物的階段結(jié)構(gòu),疾病的風(fēng)險可能被遠(yuǎn)遠(yuǎn)低估了甚至能得到滅絕和持久這樣相反的結(jié)論,也就是說有無階段結(jié)構(gòu)對于布魯氏菌病流行與否以及流行的規(guī)模具有顯著的影響.

圖5 (a) 無階段結(jié)構(gòu)模型各倉室的時間序列圖;(b) 階段結(jié)構(gòu)模型各倉室的時間序列圖;(c) 有無階段結(jié)構(gòu)I(t)的時間序列圖,所有參數(shù)值保持一致,無階段結(jié)構(gòu)的初始條件為: S(θ)=40,H(θ)=6;I(θ)=5,B(θ)=10,階段結(jié)構(gòu)的初始條件為:S1(θ)=20,S2(θ)=20,H(θ)=6;I(θ)=5,B(θ)=10,其中θ ∈[?0.2,0]