具有媒體報(bào)道和時(shí)滯效應(yīng)的SEIR模型動(dòng)力學(xué)分析

劉志華， 曹 慧

(陜西科技大學(xué) 數(shù)學(xué)與數(shù)據(jù)科學(xué)學(xué)院, 陜西 西安 710021)

0 引言

隨著信息技術(shù)的快速發(fā)展和進(jìn)步,媒體報(bào)道對控制疾病傳播起到了積極的作用,以2020年全球暴發(fā)的新冠疫情為例,媒體就疫情數(shù)據(jù)的實(shí)時(shí)報(bào)道引發(fā)了人們自覺關(guān)注疫情動(dòng)態(tài)并做出相應(yīng)的防護(hù)措施來保護(hù)自身,如勤洗手、少出門、遠(yuǎn)離密集人群等等.

近幾年來,已有不少學(xué)者借助動(dòng)力學(xué)模型研究媒體報(bào)道對疾病傳播的影響[1-6].本文將討論媒體報(bào)道對肺結(jié)核傳播的影響.

事實(shí)上,肺結(jié)核存在明顯的潛伏期,并且在潛伏期內(nèi),感染者沒有不具有傳染性,也不出現(xiàn)癥狀,很難被發(fā)現(xiàn).一旦感染者因某種原因?qū)е旅庖吡ο陆稻蜁?huì)發(fā)病成為肺結(jié)核病人.肺結(jié)核病人有癥狀和傳染性,易于發(fā)現(xiàn).而真正影響人們的是被媒體所報(bào)道出來的出現(xiàn)癥狀的發(fā)病人數(shù).因此,給出以下模型

(1)

其中:Λ是人口的常數(shù)輸入率,β是發(fā)生率,d是自然死亡率,μ是因病死亡率,α表示潛伏感染者的發(fā)病成為肺結(jié)核病人的發(fā)病率,γ表示肺結(jié)核病人的恢復(fù)率.e-εE(t-τ)表示t-τ時(shí)刻的肺結(jié)核發(fā)病人數(shù)在t時(shí)刻被媒體報(bào)道后導(dǎo)致疾病傳染率下降的作用因子,所有參數(shù)均為非負(fù)常數(shù).令N(t)表示t時(shí)刻的總?cè)丝跀?shù),即,N(t)=S(t)+E(t)+I(t)+R(t).

1 模型分析和基本再生數(shù)

1.1 模型的適定性分析

本小節(jié)將先分析系統(tǒng)(1)解的存在唯一性和非負(fù)有界性,并給出其基本再生數(shù).

定理1給定初值φ∈C+,則對所有的t>0,模型(1)存在唯一的且非負(fù)有界的解Xt.

證明：對于任意初值φ=(φ1,φ2,φ3,φ4)∈C+,定義

f(t,φ)=

其中f(t,φ)在(t,φ)∈R+×C+中是連續(xù)的,且函數(shù)f(t,φ)在C+的每個(gè)有界子集滿足Lipschitz條件,根據(jù)Hale給出的定理[7],可知模型(1)在區(qū)間t∈[0,σ],σ∈R+上存在唯一解.

假設(shè)φi(0)=0,則fi(t,φ)≥0,依據(jù)文獻(xiàn)[8]里的定理5.2.1和注釋5.2.1,可知對于t∈[0,σ],解Xt具有非負(fù)性.

將模型(1)中的四個(gè)方程相加,可得

(S(t)+E(t)+I(t)+R(t))′=

Λ-dN(t)-μI(t)≤Λ-dN(t).

所以,有

這意味著

Γ=

是模型(1)的正不變集.

證畢.

由于變量R沒有出現(xiàn)在模型(1)的前三個(gè)方程里,因此可將模型(1)簡化

(2)

顯然,模型(2)與模型(1)有相同的動(dòng)力學(xué)性態(tài).在下面的分析中,將主要討論模型(2).相應(yīng)地,模型(2)的正向不變集為

其中Lambert W(·)是一個(gè)w→wew的復(fù)值函數(shù)[9].

1.2 基本再生數(shù)

利用再生矩陣的方法[10]可得到模型(2)的基本再生數(shù)

2 平衡點(diǎn)的穩(wěn)定性

2.1 無病平衡點(diǎn)的穩(wěn)定性

定理2如果R0<1,那么模型(2)的無病平衡點(diǎn)P0是全局漸近穩(wěn)定的；如果R0>1,那么P0是不穩(wěn)定的.

證明:模型(2)在P0處的Jacobia矩陣為

可見,λ1=-d<0是一個(gè)特征根,另外兩個(gè)特征根滿足方程

λ2+(d+α+d+μ+γ)λ+(d+α)(d+μ+γ)(1-R0)=0.

這表明當(dāng)R0<1時(shí),P0是局部漸近穩(wěn)定的.反之,R0>1時(shí),P0是不穩(wěn)定的.

下面,證明無病平衡點(diǎn)的全局穩(wěn)定性.令Lyapunov函數(shù)V1(t)=αE(t)+(d+α)I(t),則V1(t)沿著模型(2)的解的全導(dǎo)數(shù)為

V1′(t)|(2)=αβe-εE(t-τ)SI-(d+α)(d+μ+γ)I≤

αβSI-(d+α)(d+μ+γ)I≤

αβS0I-(d+α)(d+μ+γ)I=

(d+α)(d+μ+γ)(R0-1)I.

證畢.

2.2 地方病平衡點(diǎn)的穩(wěn)定性

定理3如果R0>1且τ=0,那么模型(2)的地方病平衡點(diǎn)P*是全局漸近穩(wěn)定的；如果R0<1,則P*是不穩(wěn)定的.

證明:當(dāng)τ=0時(shí),模型(2)在P*處的Jacobia矩陣為

J(P*)=

其中

a=-ε(d+α)E*-(d+α),

L=(d+μ+γ)(d+α).

假設(shè)λi,i=1,2,3是矩陣J(P*)的特征根,滿足Reλ1≤Reλ2≤Reλ3.注意到d-Λ/S*<0.進(jìn)而可得detJ(P*)=(d-Λ/S*-dεE*)(d+α)(d+μ+γ)<0,所以λ1λ2λ3<0.自然對于λi,i=1,2,3的大小僅有兩種情況：Ⅰ)Reλi<0,i=1,2,3；Ⅱ) Reλ1<0≤Reλ2≤Reλ3,易看出trJ(P*)<0,故可知λ1+λ2+λ3<0,由此得出Re(λ1+λ2)<0和Re(λ1+λ3)<0.

依據(jù)文獻(xiàn)[11],J(P*)的二階加性復(fù)合矩陣為

J[2](P*)=

直接計(jì)算可得

detJ[2](P*)=

其中

K=Λ/S*+ε(d+α)E*+(d+α)>0,L>0.

根據(jù)二階加性矩陣的性質(zhì),J[2](P*)的特征根為λi+λj,1≤i1時(shí),P*是局部漸近穩(wěn)定的.

下面,證明地方病平衡點(diǎn)的全局穩(wěn)定性.令Lyapunov函數(shù)

則V2(t)沿著模型(2)的全導(dǎo)數(shù)為

V2′(t)|(2)=

證畢.

3 Hopf分支的存在性

本小節(jié)將討論當(dāng)τ>0時(shí),模型(2)在地方病平衡點(diǎn)處發(fā)生Hopf分支的充分條件.

3.1 Hopf分支的存在性

當(dāng)τ>0時(shí),模型(2)在平衡點(diǎn)P*處的Jacobia矩陣為

J(P*)=

其特征方程為

λ3+A2λ2+A1λ+A0+(B2λ2+B1λ+B0)e-λ τ=0.

(3)

其中

B2=ε(d+α)E*,

B1=ε(d+α)(d+d+μ+γ)E*,

B0=dε(d+α)(d+μ+γ)E*.

且Ai>0,Bi>0,i=0,1,2.令λ=iω(ω>0)為方程(3)的純虛根,代入可得

(4)

因此

進(jìn)一步有

ω6+p3ω4+p2ω2+p1=0,

(5)

令x=ω2,則方程(5)被重新簡化為

F(x)=x3+p3x2+p2x+p1=0.

(6)

(7)

下面,給出Reλ(t)對τ求導(dǎo)的sign大小.

定理4

證明:將方程(3)的左右兩邊分別對λ和τ進(jìn)行求導(dǎo)

因此

(8)

另外,λ=iω(ω>0)是方程(3)的純虛根,那么iω必滿足如下等式

|(iω)3+A2(iω)2+A1(iω)+A0|2=

|B2(iω)2+B1(iω)+B0|2

即

(9)

因此,利用等式(9),(8)式可簡化為

證畢.

3.2 實(shí)根的分布

本小節(jié)依據(jù)文獻(xiàn)[3]中附錄給出的一元三次方程證明,將直接給出方程(6)實(shí)根的存在條件.

F(x)=x3+p3x2+p2x+p1=0

則對于一元三次方程F(x),令

A=p32-3p2,

B=p3p2-p2p1,

C=p22-3p3p1,

Δ=B2-4AC.

引理1

(1)當(dāng)Δ>0,且p1<0時(shí),方程F(x)有且僅有一個(gè)正實(shí)根x,滿足F′(x)>0.

(2)當(dāng)Δ=0,且p3<0,p2=0,p1=0或p1<0時(shí),方程F(x)有且僅有一個(gè)正實(shí)根x,有F′(x)>0.

(3)當(dāng)Δ<0,且p2≤0或p3>0,p2>0,p1<0或p3<0,p2>0,p1≥0時(shí),方程F(x)有且僅有一個(gè)正實(shí)根x,有F′(x)>0.

(4)當(dāng)Δ<0,且p3<0,p2>0,p1<0時(shí),方程F(x)有兩個(gè)正實(shí)根xi,有F′(xi)>0,i=1,2.

(5)當(dāng)A=B=0時(shí),方程F(x)有三個(gè)重實(shí)根.

定理5如果R0>1,且

(1)滿足引理1的(1),(2),(3)中任意一個(gè)條件,則模型(2)的地方病平衡點(diǎn)P*在τ∈[0,τ0)區(qū)間內(nèi)是局部漸近穩(wěn)定的；當(dāng)τ>τ0時(shí),P*是不穩(wěn)定的,并且模型(2)會(huì)在P*附近當(dāng)τ=τ0時(shí)經(jīng)歷Hopf分支.

4 數(shù)值模擬

本小節(jié)將給出參數(shù)值,進(jìn)行數(shù)值模擬,來補(bǔ)充證實(shí)理論的正確性.

對于模型(2),選取參數(shù)Λ=0.2,β=1,ε=6,d=0.1,μ=0.1,γ=0.1.則有R0=5.925 9>1.令初始值S0=0.5,E0=0.15,I0=0.1,S0+E0+I0≤Λ/d=2,利用Matlab計(jì)算可得地方病平衡點(diǎn)P*=(S*,E*,I*)=(0.767 5,0.136 9,0.365 2).

進(jìn)一步有,Δ=-0.007 5<0,p2=-0.083 5<0.由引理1知,方程(6)存在一個(gè)正根ω1=0.242 3,再通過(7)式直接算得τ0=11.666 6,τ1=31.117 4.

給定初始值S0，E0和I0,當(dāng)τ=8<τ0時(shí)，如圖1中(a)、(b)圖所示,模型(2)的平衡點(diǎn)P*是漸近穩(wěn)定的.當(dāng)τ=20>τ0時(shí),如圖2(a)、(b)圖所示,模型(2)在平衡點(diǎn)P*附近發(fā)生Hopf分支,出現(xiàn)周期解.

隨著τ的增加,模型(2)的平衡點(diǎn)由穩(wěn)定變?yōu)椴环€(wěn)定,呈周期性變換,如圖3所示.

(a)時(shí)間序列圖

(b)相圖圖1 取τ=8<τ0時(shí),模型(2)在平衡點(diǎn)P*處的時(shí)間序列圖和相圖

(a)時(shí)間序列圖

(b)相圖圖2 取τ=20>τ0時(shí),模型(2)在平衡點(diǎn)P*處的時(shí)間序列圖和相圖

圖3 隨著τ的增加,描述模型(2)動(dòng)力學(xué)的分岔圖

5 結(jié)論

通過對SEIR模型的討論,文章證明了當(dāng)R0<1時(shí),模型(2)的無病平衡點(diǎn)P0是全局漸近穩(wěn)定的;當(dāng)R0>1,τ=0時(shí),模型(2)的地方病平衡點(diǎn)P*是全局漸近穩(wěn)定的;當(dāng)R0>1,τ>0時(shí),P*不再穩(wěn)定；進(jìn)而,在滿足引理1的條件下,當(dāng)τ=τi,i=1,2,3…時(shí),模型(2)會(huì)發(fā)生Hopf分支.

可見,媒體報(bào)道的作用因子e-εE(t-t)引起了模型(2)發(fā)生復(fù)雜的動(dòng)力學(xué)性態(tài),即,產(chǎn)生Hopf分支.