楊洛驛, 劉玲伶
(西南石油大學(xué) 理學(xué)院, 四川 成都 610500)
Covid-19是一種新型冠狀病毒感染,在病毒大流行期間,各國(guó)都采取了一些系列的有效措施來(lái)切斷其傳播,如戴口罩、病毒傳播鏈溯源、修建方艙醫(yī)院等外部管控措施,但從長(zhǎng)遠(yuǎn)來(lái)看最有效的應(yīng)該是人類自身免疫系統(tǒng)去戰(zhàn)勝病毒,其中通過(guò)疫苗接種能夠使人群對(duì)病毒產(chǎn)生免疫力和抵抗力,減少感染率和死亡率,起到間接保護(hù)效果,從而降低疫情防控造成的經(jīng)濟(jì)資源成本.疫苗接種不僅可以保護(hù)接種者自身,還可以減少傳染源和傳播途徑,同時(shí)當(dāng)社會(huì)中的大部分人群在接種有效的疫苗后,基本可以實(shí)現(xiàn)疫情的清零,因此疫苗接種為疫情防控提供更加全面的保障[1].許多學(xué)者希望通過(guò)數(shù)學(xué)建模來(lái)探討生物模型的動(dòng)態(tài)發(fā)展[2-3],研究病毒傳播模式可以模擬當(dāng)前流行病的發(fā)展,并為制定疫情管理措施提供指導(dǎo)[4].目前科學(xué)界對(duì)其流行病學(xué)、特征、臨床表現(xiàn)等方面有了更深入的研究,由于新型冠狀病毒感染的潛伏期不同,早期使用感染-再感染(SIR)模型及其擴(kuò)展的研究可能不準(zhǔn)確,而暴露-感染-再覆蓋(SEIR)模型及其延伸更為合適[5].
SEIR模型是一種廣泛應(yīng)用于現(xiàn)代流行病學(xué)的數(shù)學(xué)模型,該模型將人群分為四個(gè)相互轉(zhuǎn)化的種群:易感者、暴露者、感染者和康復(fù)者,通過(guò)對(duì)人群狀態(tài)的轉(zhuǎn)化進(jìn)行建模分析此類模型,可以有效預(yù)測(cè)疾病在人群中的發(fā)展趨勢(shì)[6].由于疫苗接種對(duì)于傳染性疾病具有重要意義,針對(duì)之前的SEIR模型沒(méi)有考慮疫苗接種對(duì)于病毒傳播的影響,2020年Nazarimehr等[7]提出了一類帶有疫苗接種的新型冠狀病毒感染SEIR模型,改進(jìn)模型如下:
(1)
其中μ是各類人群的死亡率,β是病毒從易感人群轉(zhuǎn)換到暴露人群的速率,v是疫苗接種的恢復(fù)率,σ是由暴露人群轉(zhuǎn)移到感染人群的速率,γ從感染人群轉(zhuǎn)移到恢復(fù)人群的速率.這一類模型可以用來(lái)分析新型冠狀病毒感染在人群中的傳播和疫苗接種對(duì)疫情的影響,通過(guò)這個(gè)模型可以知道一部分易感人群可以直接通過(guò)疫苗接種轉(zhuǎn)化成為康復(fù)者,以此縮小疫情的進(jìn)一步發(fā)展[8].此類模型除了可以用作新型冠狀病毒感染的研究外,還可以應(yīng)用于其他疾病的控制和管理,如流感、麻疹、水痘等,這些疾病也是全球公共衛(wèi)生的重要問(wèn)題,疫苗接種是控制和預(yù)防這些疾病的重要手段,因此,帶有疫苗接種的SEIR模型具有廣泛的應(yīng)用前景和實(shí)際意義.
本文討論了此類帶有疫苗接種的新型冠狀病毒感染SEIR模型,首先給出了系統(tǒng)平衡點(diǎn)存在的參數(shù)條件,通過(guò)計(jì)算基本再生數(shù)給出了邊界平衡點(diǎn)和內(nèi)部平衡點(diǎn)的局部穩(wěn)定性,并進(jìn)一步構(gòu)造Lyapunov函數(shù)和變分矩陣的方法給出了系統(tǒng)的全局穩(wěn)定性.當(dāng)基本再生數(shù)跨過(guò)1時(shí),討論了系統(tǒng)邊界平衡點(diǎn)附近將發(fā)生跨臨界分岔的情形,最后結(jié)合數(shù)值模擬展示系統(tǒng)穩(wěn)定性的情況,從生物學(xué)上給出模型動(dòng)力學(xué)行為的分析.
下面討論系統(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)行簡(jiǎ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)
通過(guò)求解上述方程可得:
其中
下面計(jì)算系統(tǒng)(3)的基本再生數(shù).
在邊界平衡點(diǎn)P0處,Fi(x),Vi(x)的雅可比矩陣為:
通過(guò)計(jì)算再生矩陣R0=ρ(FV-1),得到基本再生數(shù):
(5)
下面對(duì)系統(tǒng)(3)進(jìn)行平衡點(diǎn)穩(wěn)定性分析.
定理1當(dāng)R0<1時(shí),邊界平衡點(diǎn)P0是漸近穩(wěn)定的,當(dāng)R0>1時(shí),邊界平衡點(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),此時(shí)存在一個(gè)特征根為λ1=-μ-v,分析后面關(guān)于λ的一元二次方程,化簡(jiǎn)可得:
通過(guò)對(duì)應(yīng)系數(shù)相等可得:
a=(2μ+γ+σ),
(1-R0)(μ+σ)(μ+γ).
由一元二次方程的韋達(dá)定理可知,λ2λ3=b,λ2+λ3=-a,特征多項(xiàng)式的兩個(gè)特征根λ2,λ3的具體情況如下:(1)當(dāng)R0<1時(shí),特征根λ2<0,λ3<0,邊界平衡點(diǎn)P0是漸近穩(wěn)定的;(2)當(dāng)R0=1時(shí),特征根λ2=0,λ3<0,邊界平衡點(diǎn)P0是非雙曲平衡點(diǎn);(3)當(dāng)R0>1時(shí),特征根λ2>0,λ3<0,邊界平衡點(diǎn)P0是不穩(wěn)定的,定理證畢.
下面證明內(nèi)部平衡點(diǎn)P*的穩(wěn)定性.
系統(tǒng)(3)在內(nèi)部平衡點(diǎn)P*處的雅可比矩陣為:
對(duì)應(yīng)的特征多項(xiàng)式:F(λ)=a0λ3+a1λ2+a2λ+a3,對(duì)應(yīng)系數(shù)如下:
a0=1,a1=γ+3μ+σ+v+(μ+v)(R0-1),a2=(μ+v)(γ+2μ+σ)R0,a3=(μ+v)(μ+γ)(μ+σ)(R0-1).
當(dāng)R0>1時(shí),滿足ai>0,計(jì)算可得a1a2-a3>0.因此當(dāng)R0>1時(shí),有主子式
Δ3:=a3·Δ2>0,
根據(jù)Routh-Hurwitz定理可得,內(nèi)部平衡點(diǎn)P*是漸近穩(wěn)定的,因此定理1得證.
下面進(jìn)行系統(tǒng)(3)的全局漸近穩(wěn)定性分析.
定理2當(dāng)R0<1時(shí),邊界平衡點(diǎn)P0是全局漸近穩(wěn)定的;當(dāng)R0>1時(shí),內(nèi)部平衡點(diǎn)P*是全局漸近穩(wěn)定的.
證在可行域T內(nèi),構(gòu)造Lyapunov函數(shù)
V=σ·Y+(μ+σ)·Z≥0,
當(dāng)R0<1時(shí),
dV/dt=σ·(βXZ-(μ+σ)Y)+(μ+σ)·
(σY-(μ+γ)Z)=
Z·(μ+σ)(μ+v)(R0-1)<0,
其中,當(dāng)Z=0時(shí),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*處的雅可比矩陣如下:
則對(duì)應(yīng)的第二變分矩陣如下:
其中
R3上的范數(shù)
|(μ,v,w)|=max{|μ|,|v|+|w|},
其中
η(B)≤max{g1,g2},
g1=η1(B11)+|B12|,g2=η1(B22)+|B21|,
η1是相對(duì)于L1矩陣范數(shù)的測(cè)度,而|B12|,|B21|是L1的矩陣范式.通過(guò)計(jì)算可得:
因此
下面將對(duì)系統(tǒng)(3)在R0=1處的分岔情況進(jìn)行分析.
定理3在R0=1處,系統(tǒng)(3)在邊界平衡點(diǎn)P0附近將會(huì)發(fā)生跨臨界分岔.
證將平衡點(diǎn)在P0處移到原點(diǎn),進(jìn)行坐標(biāo)平移,令
系統(tǒng)(3)轉(zhuǎn)化為如下等價(jià)系統(tǒng):
(6)
令分岔參數(shù)ξ=βμσ-(σ+μ)(γ+μ)(v+μ),可得
當(dāng)R0=1時(shí),ξ=0,雅可比矩陣如下:
特征多項(xiàng)式有特征根
λ1=-μ-v,λ2=0,λ3=-2μ-γ-σ.
對(duì)于特征根λ2=0,有對(duì)應(yīng)特征向量
上述雅可比矩陣的轉(zhuǎn)置如下:
對(duì)于上述轉(zhuǎn)置后矩陣有特征根λ=0,對(duì)應(yīng)的特征向量φ=(φ1,φ2,φ3)T=(0,σ,μ+σ)T.
令f=(f1,f2,f3),則原方程可得:
通過(guò)計(jì)算得到:
因此
所以由Sotomayor定理[12]知,在ξ=0(R0=1)處,系統(tǒng)(3)將在P0附近發(fā)生跨臨界分岔,定理3得證.
定理(5)分析的分岔現(xiàn)象表明,平衡點(diǎn)的穩(wěn)定性隨著參數(shù)取值的變化而發(fā)生突變,從而導(dǎo)致病毒傳播速度不斷增加,需要采取有效措施來(lái)控制疫情的蔓延,如限制人口流動(dòng)、加強(qiáng)衛(wèi)生防護(hù)、推廣有效疫苗接種等.同時(shí)應(yīng)密切關(guān)注病毒傳播的發(fā)展趨勢(shì),阻止或延緩分岔現(xiàn)象的出現(xiàn),避免造成更大的經(jīng)濟(jì)與人力資源損失.
本節(jié)利用Matlab軟件中“ode45”對(duì)系統(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)此時(shí)僅存在一個(gè)邊界平衡點(diǎn),且由定理2知邊界平衡點(diǎn)是全局漸近穩(wěn)定的.圖1(d)顯示,隨著時(shí)間變化,軌線最終趨于邊界平衡點(diǎn).分析各類人群可得,易感人群(S)隨著時(shí)間先下降并趨于穩(wěn)定(見(jiàn)圖1(a));由圖1(b)、(c)可知暴露人群(E)和感染人群(I)隨時(shí)間不斷降低,疫情得到了有效控制.但需要注意的是,即使感染人群已經(jīng)降低至零,仍然存在一部分易感人群,因此需要加強(qiáng)疫情監(jiān)測(cè)和防控工作,及時(shí)發(fā)布疫情信息,提高公眾疫情防控意識(shí),及時(shí)采取必要的隔離和治療措施,同時(shí)應(yīng)加強(qiáng)社會(huì)宣傳和防疫教育,提高公眾對(duì)病毒傳播途徑、防護(hù)措施等方面的了解,注意個(gè)人防護(hù)和糾正不良衛(wèi)生習(xí)慣,增強(qiáng)身體免疫力,更好地保護(hù)自己和他人.
選取參數(shù)μ=0.012,β=0.560,v=0.010,σ=0.050,γ=0.003,由基本再生數(shù)(5)可得R0>1,由引理1可知此時(shí)存在一個(gè)內(nèi)部平衡點(diǎn).且由定理2可知內(nèi)部平衡點(diǎn)是全局漸近穩(wěn)定的.圖2(d)可以看出,隨時(shí)間變化,軌線趨于內(nèi)部平衡點(diǎn);觀察各類人群,由圖2(a)可知,易感人群(S)隨時(shí)間不斷減少并逐步增加至定值;圖2(b)顯示,暴露人群(E)隨時(shí)間先增加再減少至定值;圖2(c)可知感染人群(I)隨時(shí)間迅速增長(zhǎng)并緩慢減少到達(dá)飽和狀態(tài).通過(guò)分析可知,在疫情傳播初期病毒的傳播速度非常迅速,容易導(dǎo)致疫情快速擴(kuò)大.因此,需要盡快采取措施切斷病毒的傳播途徑,比如加強(qiáng)個(gè)人防護(hù)、避免密閉空間的聚集、加強(qiáng)社交距離的管理等措施.此外,針對(duì)感染者和密切接觸者,應(yīng)及時(shí)進(jìn)行隔離和檢測(cè)管理措施,同時(shí)可以看到隨著時(shí)間的推移,感染者數(shù)量逐漸達(dá)到一個(gè)飽和峰值,這個(gè)階段應(yīng)當(dāng)積極協(xié)調(diào)資源,做好感染后的治療措施,加強(qiáng)完善醫(yī)療衛(wèi)生系統(tǒng)的建設(shè),避免造成二次感染.
(a)易感人群
(b)暴露人群
(c)感染人群
(d)解
(a)易感人群的軌跡圖
(b)暴露人群
(c)感染人群
(d)解
內(nèi)江師范學(xué)院學(xué)報(bào)2023年12期