苗 慧，李淑萍

(中北大學(xué) 理學(xué)院，山西 太原 030051)

0 引 言

從2019年起，由嚴(yán)重急性呼吸綜合征冠狀病毒2(SARS-CoV-2)引起的新冠肺炎在世界各地大規(guī)模流行.據(jù)世界衛(wèi)生組織WHO報告數(shù)據(jù)來看，截至2022年1月17日,全世界累計確診329 109 769例，死亡555 3630例.全球新冠肺炎病例不斷增加，全面了解它的傳播途徑，采取對應(yīng)的防控措施非常重要.研究表明，呼吸道飛沫和直接接觸是新冠病毒傳播的主要途徑[1-2]，且新冠病毒在物體表面也能存活數(shù)天[3]，因此，人類通過間接接觸被污染的物體表面也可能會導(dǎo)致感染[4-5].

專家學(xué)者對不同地區(qū)新冠肺炎的傳播進(jìn)行了數(shù)學(xué)建模.Wang Liyan等在揚(yáng)州疫情的基礎(chǔ)上考慮了潛伏期和治愈期所服從的新密度函數(shù)，以及感染模型中的疫苗接種效應(yīng)，研究了長潛伏期新冠肺炎的傳播機(jī)制，認(rèn)為通過足夠的隔離和高接種疫苗措施可促使疫情消失[6].Kolebaje O T等提出了一種結(jié)合非藥物干預(yù)措施的新冠肺炎SEIIAQHR數(shù)學(xué)模型，其中S代表易感者，E代表潛伏者，I代表感染者，IA代表無癥狀感染者，Q代表隔離檢疫者，H代表住院者，R代表恢復(fù)者，并將此模型用于研究一些非洲國家的累計確診病例數(shù)的時間序列演變[7].Shen Zhonghua等考慮了巴基斯坦的新冠肺炎報告病例，建立了帶有疫苗接種V倉室的SVEAIR模型(A為無癥狀感染倉室)，并進(jìn)行最優(yōu)控制分析[8].Abioye A I等建立了帶有隔離倉室的SEIQR模型，并結(jié)合戴口罩、積極篩查檢測、控制社交距離等嚴(yán)格的控制策略，預(yù)測尼日利亞的新冠肺炎將被根除[9].

上述模型有的考慮了隔離檢疫，有的考慮了疫苗接種，但沒有考慮環(huán)境病毒對疫情傳播的影響.自新冠疫情爆發(fā)以來，我國政府采取了一系列隔離防控措施，如通過建立火神山醫(yī)院集中隔離收治患者[10]，并取得了很大的成效.因此，本文建立了一個隔離倉室Q的SEIQR模型.另外，近期的深圳疫情很可能是由國外入境的新冠病毒污染物品引起的[11]，本文將研究環(huán)境病毒W(wǎng)對新冠肺炎傳播的影響，以期為政府部門制定防控策略提供理論依據(jù).

1 模型的建立

S(t)為t時刻易感者數(shù)量，E(t)為t時刻潛伏者數(shù)量，I(t)為t時刻感染者數(shù)量，Q(t)為t時刻隔離檢疫者數(shù)量，R(t)為t時刻恢復(fù)者數(shù)量，W(t) 為t時刻環(huán)境中的新冠病毒數(shù)量，模型倉室圖如圖1 所示.圖1 中，Λ為總?cè)丝诘妮斎肼剩?為潛伏者對易感者的傳染率，β2為染病者對易感者的傳染率，β3為環(huán)境中的新冠病毒對易感者的傳染率，μ為自然死亡率，k為潛伏者到染病者的轉(zhuǎn)化率，b為染病者隔離檢疫的比例系數(shù)，α為染病者的因病死亡率，γ為染病者的恢復(fù)率，δ為隔離后的康復(fù)率，m為潛伏者向外界環(huán)境釋放病毒的速率，μ1為環(huán)境中新冠病毒的自然死亡率.模型為

圖1 SEIQRW流程圖Fig.1 Flow chart of SEIQRW transmission

(1)

由N=S+E+I+Q+R得

(2)

Ω={(S,E,I,Q,R,W)|S,E,I,Q,R,W≥0,

(3)

2 平衡點(diǎn)和基本再生數(shù)的計算

利用下一代矩陣法[12]，系統(tǒng)(1)中受感染的相關(guān)倉室可寫為

式中：X=(E,I,Q,W)T；F為新感染的矩陣；V為倉室個體轉(zhuǎn)移矩陣，則

在無病平衡點(diǎn)E0處對(E,I,Q,W)求偏導(dǎo)，可得

計算得

因此，基本再生數(shù)為

為求解正平衡點(diǎn)，令模型(1)右端為0，則

(4)

當(dāng)R0>1時，正平衡點(diǎn)為E1=(S*,E*,I*,Q*,R*,W*).其中，

3 平衡點(diǎn)的穩(wěn)定性分析

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

定理 1當(dāng)R0<1時，系統(tǒng)(1)的無病平衡點(diǎn)E0在Ω內(nèi)是局部漸近穩(wěn)定的；當(dāng)R0>1時,無病平衡點(diǎn)不穩(wěn)定.

證明系統(tǒng)(1)在E0的Jacobian矩陣為

J|E0的特征方程為

(λ+μ)2(λ+μ+α+δ)(λ3+a1λ2+

a2λ+a3)=0.

(5)

顯而易見，λ1,λ2=-μ,λ3=-(μ+α+δ).

其中，

μ1+b+μ+α+γ>0,

a2=(μ+k)(b+μ+α+γ)+μ1(μ+k)+

μ1(b+μ+α+γ)>0,

a3=μ1(μ+k)(b+μ+α+γ)(1-R0)>0,

μ1(b+μ+α+γ)2>0.

由Hurwitz判據(jù)[13]可知，無病平衡點(diǎn)E0在可行域內(nèi)局部漸近穩(wěn)定.

定理 2當(dāng)R0<1時，系統(tǒng)(1)的無病平衡點(diǎn)E0在Ω內(nèi)是全局漸近穩(wěn)定的.

證明可以將原系統(tǒng)(1)的染病倉室重新表示為

x′=(F-V)x-f(x,y),

(6)

式中：x=(E,I,Q,W)T為染病倉室中的人群；y=(S,R)T為非染病倉室中的人群，f(x,y)=(0,0,0,0)T.

運(yùn)用矩陣?yán)碚撝械腜erron特征向量[14]，非負(fù)矩陣V-1F的左特征向量ωT對應(yīng)的特征值為ρ(V-1F)=ρ(FV-1)=R0，可求得ωT=(β1,β2,0,β3).若f(x,y)≥0,F≥0,V-1≥0，則定義Lyapunov函數(shù)為

Q=ωTV-1x=

當(dāng)R0<1時，

Q′=ωTV-1x′=

(R0-1)ωTx-ωTV-1f(x,y)=

(R0-1)(β1E+β2I+β3W)≤0.

所以，當(dāng)R0<1時，無病平衡點(diǎn)X0在Ω內(nèi)全局漸近穩(wěn)定.

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

b1=A+B+C+μ+μ1,

b2=A(B+C)+(A+C)μ1+(B+C)μ+μμ1+

b3=ABC+(AB+AC)μ1+Cμμ1+Aμ1β1S*+

b4=ABCμ1+ACμ1β1S*,

A=β1E*+β2I*+β3W*,

B=μ+k-β1S*=[kμ1β2(μ+k)+

mβ3C(μ+k)]/Ζ,

C=b+μ+α+γ,

Z=μ1β1C+kμ1β2+mβ3C.

證明系統(tǒng)(1)在E1處的Jacobian矩陣為

J|E1的特征方程為

(λ+μ)(λ+μ+α+δ)(λ4+b1λ3+b2λ2+

b3λ+b4)=0.

(7)

顯然，λ1=-μ,λ2=-(μ+α+δ).b1,b2,b3,b4>0，且有

b1b2-b3=(A+B)a2+ABμ+ACμ+AC2+

μ1(A+C)(C+μ+μ1)+μ(B+C)(C+μ+

μ1)+μμ1(μ+μ1)+Aμβ1S*+(C+μ1)·

4 最優(yōu)控制策略

為找到合理有效的控制策略，本文采用最優(yōu)控制理論[15]，并運(yùn)用龐特里亞金極大值原理[16]找到模型中最優(yōu)控制所需要的條件.

在系統(tǒng)(1)中引入3個控制項u1,u2,u3，將控制系統(tǒng)表示為

(8)

式中：u1為減少易感者與新冠病毒的直接接觸和間接接觸的隔離策略；u2為增強(qiáng)感染者的治療恢復(fù)率的治療策略；u3為消滅環(huán)境中病毒的環(huán)境衛(wèi)生策略.

控制函數(shù)集為

U={U(·)∈L1([0,T]；R3)|0

umax<1,?t∈[0,T]}.

為考慮感染者的治療成本和控制成本，將目標(biāo)函數(shù)定義為

(9)

龐特里亞金極大值原理是將式(8)和(9)轉(zhuǎn)化為一個關(guān)于控制變量的最小化的哈密頓問題，構(gòu)造哈密頓函數(shù)，即

(10)

式中：λi(i=1,2,…,6)為伴隨變量.

(11)

那么，最優(yōu)控制為

證明由伴隨系統(tǒng)(11)得

根據(jù)最優(yōu)化條件

得最優(yōu)控制值為

5 數(shù)值模擬

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

假設(shè)Λ=100,β1=0.000 2,β2=0.000 1,β3=0.000 05,μ=0.069,μ1=0.06,γ=0.04,k=0.69,b=0.7,α=0.006 6,m=0.009,δ=0.83，可得R0=0.56<1，如圖2(a)所示，驗證了無病平衡點(diǎn)是全局漸近穩(wěn)定的.

(a) R0<1

假設(shè)β1=0.000 9,β2=0.000 2,β3=0.000 05,m=0.09,b=0.07,δ=0.63，可得R0=3.28>1，如圖2(b)所示，故地方病平衡點(diǎn)是全局漸近穩(wěn)定的.

5.2 數(shù)據(jù)擬合

基于政府在新冠肺炎疫情前期采取隔離防控措施時感染者的數(shù)據(jù)以及部分參數(shù)估計[17]，得到k=1/5.2,γ=1/7,α=0.043,δ=1/10,并且假設(shè)Λ=3 000,β1=0.65×10-5,β2=0.85×10-5,β3=0.59×10-6,μ=0.06,μ1=0.056,m=0.89,b=0.07.運(yùn)用Matlab模擬對比了有無考慮環(huán)境病毒時估計的染病者人數(shù)曲線，如圖3 所示，發(fā)現(xiàn)帶病毒項的系統(tǒng)與實(shí)際感染人數(shù)更相符.

圖3 帶病毒、不帶病毒的模型與患病人數(shù)擬合Fig.3 Fitting model with and without virus to number of sick people

5.3 敏感性分析

運(yùn)用PRCC法進(jìn)行敏感性分析，可以更好地觀察各個參數(shù)與基本再生數(shù)R0的相關(guān)性，如圖4 所示.可知參數(shù)Λ,β1,β2,β3,m與R0呈正相關(guān)，μ,μ1,k,b,α,γ與R0呈負(fù)相關(guān).這也說明了直接和間接接觸都加速了病毒傳播，而控制接觸和消滅環(huán)境病毒可降低基本再生數(shù)，從而減少感染人數(shù).

圖4 R0與參數(shù)的相關(guān)性Fig.4 Correlation between R0 and parameters

5.4 最優(yōu)控制

利用前后掃描法和四階龍格庫塔公式進(jìn)行模擬，為控制成本，節(jié)約醫(yī)療資源，令u1=0.8,u2=0.7,u3=0.6，得到最優(yōu)控制函數(shù)圖，如圖5(a) 所示.令Λ=100,β1=0.000 2,β2=0.000 1,β3=0.000 05,μ=0.069,μ1=0.06,m=0.009,k=0.69,b=0.7,α=0.006 6,γ=0.04,δ=0.83，計算得R0=0.56<1.如圖5(b)和(c)，當(dāng)加入控制時，感染人數(shù)和環(huán)境中病毒數(shù)量在較短時間內(nèi)迅速減少為0，由此可知，最優(yōu)控制策略可以阻斷新冠肺炎的傳播.

(a) 控制變量u(t)與t的函數(shù)關(guān)系

6 結(jié) 論

本文建立了一個包括隔離和環(huán)境中新冠病毒倉室的模型，通過動力學(xué)分析，得到基本再生數(shù)R0、無病平衡點(diǎn)和正平衡點(diǎn)的穩(wěn)定性及最優(yōu)控制策略.基于數(shù)值模擬對理論分析進(jìn)行了驗證，參數(shù)的敏感性分析和數(shù)據(jù)擬合進(jìn)一步證實(shí)隔離是疫情防控必不可少的措施.將環(huán)境病毒考慮在內(nèi)可使得模型曲線與患病人數(shù)更加接近，強(qiáng)調(diào)了環(huán)境病毒在疫情傳播過程是不可忽視的.如果不在模型中納入環(huán)境間接傳播，會導(dǎo)致低估疫情的嚴(yán)重程度.研究結(jié)果為世衛(wèi)組織建議的洗手和表面清潔和消毒的做法提供了強(qiáng)有力的支撐，也為我國采取及時隔離的政策提供了理論依據(jù).