楊洛驛, 劉玲伶

(西南石油大學(xué) 理學(xué)院, 四川 成都 610500)

0 引言

Covid-19是一種新型冠狀病毒感染,在病毒大流行期間,各國都采取了一些系列的有效措施來切斷其傳播,如戴口罩、病毒傳播鏈溯源、修建方艙醫(yī)院等外部管控措施,但從長遠(yuǎn)來看最有效的應(yīng)該是人類自身免疫系統(tǒng)去戰(zhàn)勝病毒,其中通過疫苗接種能夠使人群對病毒產(chǎn)生免疫力和抵抗力,減少感染率和死亡率,起到間接保護(hù)效果,從而降低疫情防控造成的經(jīng)濟(jì)資源成本.疫苗接種不僅可以保護(hù)接種者自身,還可以減少傳染源和傳播途徑,同時當(dāng)社會中的大部分人群在接種有效的疫苗后,基本可以實(shí)現(xiàn)疫情的清零,因此疫苗接種為疫情防控提供更加全面的保障[1].許多學(xué)者希望通過數(shù)學(xué)建模來探討生物模型的動態(tài)發(fā)展[2-3],研究病毒傳播模式可以模擬當(dāng)前流行病的發(fā)展,并為制定疫情管理措施提供指導(dǎo)[4].目前科學(xué)界對其流行病學(xué)、特征、臨床表現(xiàn)等方面有了更深入的研究,由于新型冠狀病毒感染的潛伏期不同,早期使用感染-再感染(SIR)模型及其擴(kuò)展的研究可能不準(zhǔn)確,而暴露-感染-再覆蓋(SEIR)模型及其延伸更為合適[5].

SEIR模型是一種廣泛應(yīng)用于現(xiàn)代流行病學(xué)的數(shù)學(xué)模型,該模型將人群分為四個相互轉(zhuǎn)化的種群:易感者、暴露者、感染者和康復(fù)者,通過對人群狀態(tài)的轉(zhuǎn)化進(jìn)行建模分析此類模型,可以有效預(yù)測疾病在人群中的發(fā)展趨勢[6].由于疫苗接種對于傳染性疾病具有重要意義,針對之前的SEIR模型沒有考慮疫苗接種對于病毒傳播的影響,2020年Nazarimehr等[7]提出了一類帶有疫苗接種的新型冠狀病毒感染SEIR模型,改進(jìn)模型如下:

(1)

其中μ是各類人群的死亡率,β是病毒從易感人群轉(zhuǎn)換到暴露人群的速率,v是疫苗接種的恢復(fù)率,σ是由暴露人群轉(zhuǎn)移到感染人群的速率,γ從感染人群轉(zhuǎn)移到恢復(fù)人群的速率.這一類模型可以用來分析新型冠狀病毒感染在人群中的傳播和疫苗接種對疫情的影響,通過這個模型可以知道一部分易感人群可以直接通過疫苗接種轉(zhuǎn)化成為康復(fù)者,以此縮小疫情的進(jìn)一步發(fā)展[8].此類模型除了可以用作新型冠狀病毒感染的研究外,還可以應(yīng)用于其他疾病的控制和管理,如流感、麻疹、水痘等,這些疾病也是全球公共衛(wèi)生的重要問題,疫苗接種是控制和預(yù)防這些疾病的重要手段,因此,帶有疫苗接種的SEIR模型具有廣泛的應(yīng)用前景和實(shí)際意義.

本文討論了此類帶有疫苗接種的新型冠狀病毒感染SEIR模型,首先給出了系統(tǒng)平衡點(diǎn)存在的參數(shù)條件,通過計(jì)算基本再生數(shù)給出了邊界平衡點(diǎn)和內(nèi)部平衡點(diǎn)的局部穩(wěn)定性,并進(jìn)一步構(gòu)造Lyapunov函數(shù)和變分矩陣的方法給出了系統(tǒng)的全局穩(wěn)定性.當(dāng)基本再生數(shù)跨過1時,討論了系統(tǒng)邊界平衡點(diǎn)附近將發(fā)生跨臨界分岔的情形,最后結(jié)合數(shù)值模擬展示系統(tǒng)穩(wěn)定性的情況,從生物學(xué)上給出模型動力學(xué)行為的分析.

1 平衡點(diǎn)存在性和局部穩(wěn)定性分析

下面討論系統(tǒng)(3)平衡點(diǎn)的存在性和穩(wěn)定性.

考慮生物意義可知各參數(shù)均大于0,其中系統(tǒng)(1)人數(shù)N(t)=S(t)+E(t)+I(t)+R(t),且dN/dt≡0,可知總?cè)藬?shù)為定值,可將系統(tǒng)(1)進(jìn)行簡化:令變換

X=S/N,Y=E/N,Z=I/N,H=R/N,

則將系統(tǒng)(1)轉(zhuǎn)化為:

(2)

由系統(tǒng)(1)的人群關(guān)系可得

X(t)+Y(t)+Z(t)+H(t)=1,

H(t)=1-X(t)-Y(t)-Z(t).

利用人口總數(shù)是常數(shù)的假設(shè),則可將上述系統(tǒng)(2)轉(zhuǎn)化為如下三維系統(tǒng):

(3)

該模型的可行域?yàn)?

證系統(tǒng)平衡點(diǎn)滿足如下條件:

(4)

通過求解上述方程可得:

其中

下面計(jì)算系統(tǒng)(3)的基本再生數(shù).

在邊界平衡點(diǎn)P0處,Fi(x),Vi(x)的雅可比矩陣為:

通過計(jì)算再生矩陣R0=ρ(FV-1),得到基本再生數(shù):

(5)

下面對系統(tǒng)(3)進(jìn)行平衡點(diǎn)穩(wěn)定性分析.

定理1當(dāng)R0<1時,邊界平衡點(diǎn)P0是漸近穩(wěn)定的,當(dāng)R0>1時,邊界平衡點(diǎn)P0是不穩(wěn)定的,內(nèi)部平衡點(diǎn)P*是漸近穩(wěn)定的.

證邊界平衡點(diǎn)P0的雅可比矩陣如下:

系統(tǒng)(3)在邊界平衡點(diǎn)P0處的特征多項(xiàng)式為(λ+μ+v)(λ2+aλ+b),此時存在一個特征根為λ1=-μ-v,分析后面關(guān)于λ的一元二次方程,化簡可得:

通過對應(yīng)系數(shù)相等可得:

a=(2μ+γ+σ),

(1-R0)(μ+σ)(μ+γ).

由一元二次方程的韋達(dá)定理可知,λ2λ3=b,λ2+λ3=-a,特征多項(xiàng)式的兩個特征根λ2,λ3的具體情況如下:(1)當(dāng)R0<1時,特征根λ2<0,λ3<0,邊界平衡點(diǎn)P0是漸近穩(wěn)定的;(2)當(dāng)R0=1時,特征根λ2=0,λ3<0,邊界平衡點(diǎn)P0是非雙曲平衡點(diǎn);(3)當(dāng)R0>1時,特征根λ2>0,λ3<0,邊界平衡點(diǎn)P0是不穩(wěn)定的,定理證畢.

下面證明內(nèi)部平衡點(diǎn)P*的穩(wěn)定性.

系統(tǒng)(3)在內(nèi)部平衡點(diǎn)P*處的雅可比矩陣為:

對應(yīng)的特征多項(xiàng)式:F(λ)=a0λ3+a1λ2+a2λ+a3,對應(yīng)系數(shù)如下:

a0=1,a1=γ+3μ+σ+v+(μ+v)(R0-1),a2=(μ+v)(γ+2μ+σ)R0,a3=(μ+v)(μ+γ)(μ+σ)(R0-1).

當(dāng)R0>1時,滿足ai>0,計(jì)算可得a1a2-a3>0.因此當(dāng)R0>1時,有主子式

Δ3:=a3·Δ2>0,

根據(jù)Routh-Hurwitz定理可得,內(nèi)部平衡點(diǎn)P*是漸近穩(wěn)定的,因此定理1得證.

2 全局穩(wěn)定性

下面進(jìn)行系統(tǒng)(3)的全局漸近穩(wěn)定性分析.

定理2當(dāng)R0<1時,邊界平衡點(diǎn)P0是全局漸近穩(wěn)定的;當(dāng)R0>1時,內(nèi)部平衡點(diǎn)P*是全局漸近穩(wěn)定的.

證在可行域T內(nèi),構(gòu)造Lyapunov函數(shù)

V=σ·Y+(μ+σ)·Z≥0,

當(dāng)R0<1時,

dV/dt=σ·(βXZ-(μ+σ)Y)+(μ+σ)·

(σY-(μ+γ)Z)=

Z·(μ+σ)(μ+v)(R0-1)<0,

其中,當(dāng)Z=0時,dV/dt=0,由Lasalle不變集原理可得,邊界平衡點(diǎn)P0是全局漸近穩(wěn)定的.

下面證明內(nèi)部平衡點(diǎn)P*的全局穩(wěn)定性.

利用文獻(xiàn)[10]中變分矩陣的方法證明內(nèi)部平衡點(diǎn)P*全局漸近穩(wěn)定性,系統(tǒng)(3)內(nèi)部平衡點(diǎn)P*處的雅可比矩陣如下:

則對應(yīng)的第二變分矩陣如下:

其中

R3上的范數(shù)

|(μ,v,w)|=max{|μ|,|v|+|w|},

其中

η(B)≤max{g1,g2},

g1=η1(B11)+|B12|,g2=η1(B22)+|B21|,

η1是相對于L1矩陣范數(shù)的測度,而|B12|,|B21|是L1的矩陣范式.通過計(jì)算可得:

因此

3 跨臨界分岔

下面將對系統(tǒng)(3)在R0=1處的分岔情況進(jìn)行分析.

定理3在R0=1處,系統(tǒng)(3)在邊界平衡點(diǎn)P0附近將會發(fā)生跨臨界分岔.

證將平衡點(diǎn)在P0處移到原點(diǎn),進(jìn)行坐標(biāo)平移,令

系統(tǒng)(3)轉(zhuǎn)化為如下等價系統(tǒng):

(6)

令分岔參數(shù)ξ=βμσ-(σ+μ)(γ+μ)(v+μ),可得

當(dāng)R0=1時,ξ=0,雅可比矩陣如下:

特征多項(xiàng)式有特征根

λ1=-μ-v,λ2=0,λ3=-2μ-γ-σ.

對于特征根λ2=0,有對應(yīng)特征向量

上述雅可比矩陣的轉(zhuǎn)置如下:

對于上述轉(zhuǎn)置后矩陣有特征根λ=0,對應(yīng)的特征向量φ=(φ1,φ2,φ3)T=(0,σ,μ+σ)T.

令f=(f1,f2,f3),則原方程可得:

通過計(jì)算得到:

因此

所以由Sotomayor定理[12]知,在ξ=0(R0=1)處,系統(tǒng)(3)將在P0附近發(fā)生跨臨界分岔,定理3得證.

定理(5)分析的分岔現(xiàn)象表明,平衡點(diǎn)的穩(wěn)定性隨著參數(shù)取值的變化而發(fā)生突變,從而導(dǎo)致病毒傳播速度不斷增加,需要采取有效措施來控制疫情的蔓延,如限制人口流動、加強(qiáng)衛(wèi)生防護(hù)、推廣有效疫苗接種等.同時應(yīng)密切關(guān)注病毒傳播的發(fā)展趨勢,阻止或延緩分岔現(xiàn)象的出現(xiàn),避免造成更大的經(jīng)濟(jì)與人力資源損失.

4 數(shù)值模擬和生物意義

本節(jié)利用Matlab軟件中“ode45”對系統(tǒng)(3)進(jìn)行數(shù)值模擬.

選取參數(shù)μ=0.012,β=0.090,v=0.100,σ=0.038,γ=0.090,由基本再生數(shù)(5)計(jì)算可得R0<1,由引理1可得系統(tǒng)(3)此時僅存在一個邊界平衡點(diǎn),且由定理2知邊界平衡點(diǎn)是全局漸近穩(wěn)定的.圖1(d)顯示,隨著時間變化,軌線最終趨于邊界平衡點(diǎn).分析各類人群可得,易感人群(S)隨著時間先下降并趨于穩(wěn)定(見圖1(a));由圖1(b)、(c)可知暴露人群(E)和感染人群(I)隨時間不斷降低,疫情得到了有效控制.但需要注意的是,即使感染人群已經(jīng)降低至零,仍然存在一部分易感人群,因此需要加強(qiáng)疫情監(jiān)測和防控工作,及時發(fā)布疫情信息,提高公眾疫情防控意識,及時采取必要的隔離和治療措施,同時應(yīng)加強(qiáng)社會宣傳和防疫教育,提高公眾對病毒傳播途徑、防護(hù)措施等方面的了解,注意個人防護(hù)和糾正不良衛(wèi)生習(xí)慣,增強(qiáng)身體免疫力,更好地保護(hù)自己和他人.

選取參數(shù)μ=0.012,β=0.560,v=0.010,σ=0.050,γ=0.003,由基本再生數(shù)(5)可得R0>1,由引理1可知此時存在一個內(nèi)部平衡點(diǎn).且由定理2可知內(nèi)部平衡點(diǎn)是全局漸近穩(wěn)定的.圖2(d)可以看出,隨時間變化,軌線趨于內(nèi)部平衡點(diǎn);觀察各類人群,由圖2(a)可知,易感人群(S)隨時間不斷減少并逐步增加至定值;圖2(b)顯示,暴露人群(E)隨時間先增加再減少至定值;圖2(c)可知感染人群(I)隨時間迅速增長并緩慢減少到達(dá)飽和狀態(tài).通過分析可知,在疫情傳播初期病毒的傳播速度非常迅速,容易導(dǎo)致疫情快速擴(kuò)大.因此,需要盡快采取措施切斷病毒的傳播途徑,比如加強(qiáng)個人防護(hù)、避免密閉空間的聚集、加強(qiáng)社交距離的管理等措施.此外,針對感染者和密切接觸者,應(yīng)及時進(jìn)行隔離和檢測管理措施,同時可以看到隨著時間的推移,感染者數(shù)量逐漸達(dá)到一個飽和峰值,這個階段應(yīng)當(dāng)積極協(xié)調(diào)資源,做好感染后的治療措施,加強(qiáng)完善醫(yī)療衛(wèi)生系統(tǒng)的建設(shè),避免造成二次感染.

(a)易感人群

(b)暴露人群

(c)感染人群

(d)解

(a)易感人群的軌跡圖

(b)暴露人群

(c)感染人群

(d)解