魏慶棟，王 忠，湯家軍，陳海燕

（1.火箭軍工程大學(xué)，西安 710038；2.中國礦業(yè)大學(xué)（北京），北京 100083）

自新冠疫情防控進(jìn)入常態(tài)化以來，使用合適的數(shù)學(xué)模型來描述和預(yù)測病毒傳播情況，評(píng)估疫情防控措施效果已經(jīng)成為學(xué)界共識(shí)。本文通過研究新冠肺炎傳播的動(dòng)力學(xué)特性，引入時(shí)滯參數(shù)，結(jié)合分?jǐn)?shù)階模型，構(gòu)建了一種描述新冠肺炎傳播狀況的分?jǐn)?shù)階群體時(shí)滯模型。

1 模型的建立與分析

1.1 分?jǐn)?shù)階模型的建立

根據(jù)COVID-19的傳播特性，在傳染病動(dòng)力學(xué)SIR模型中引入潛伏期患者，即已經(jīng)感染病毒但是尚未表現(xiàn)出癥狀的患者，也可以稱為無癥狀患者，因此在模型中增加時(shí)滯參數(shù)，根據(jù)現(xiàn)有經(jīng)驗(yàn)，在經(jīng)過一段時(shí)間的潛伏后，會(huì)有一定比例的潛伏患者最終發(fā)展成為確診病例[1]，這部分人群對(duì)最終結(jié)果的影響離不開潛伏周期這一因素，用參數(shù)τ來表示，在完成整數(shù)階模型后，修改得到下列分?jǐn)?shù)階模型

上述系統(tǒng)當(dāng)t=0時(shí)，會(huì)滿足下述的初值條件

其中，S（t）為易感人群t時(shí)刻的人數(shù)，A（t）表示在t時(shí)刻無癥狀感染者的人數(shù)，I（t）為在t時(shí)刻有癥狀的感染者的人數(shù)，R（t）為恢復(fù)健康的人數(shù)，λ為正常人群中的自然增長速率下的增長人數(shù)，θ為無癥狀感染者轉(zhuǎn)換為有癥狀感染者的概率，β1、β2分別為正常人群和無癥狀感染者與有癥狀感染者接觸后被感染的概率，δA、δI分別為無癥狀被感染人群和有癥狀被感染人群的恢復(fù)率，d、dA、dI分別為正常人群、無癥狀感染者和有癥狀感染者的死亡率。在本系統(tǒng)中引入時(shí)滯參數(shù)τ來模擬模型的潛伏期參數(shù)。

1.2 模型的分析

根據(jù)系統(tǒng)1（公式（1）所示，下同），可以得出2個(gè)平衡點(diǎn)，一個(gè)無病平衡點(diǎn)E0=（x0，0，0，0），其中，一個(gè)地方病平衡點(diǎn)E*=（S*，A*，I*，R*），該平衡點(diǎn)滿足

在對(duì)上述2個(gè)平衡點(diǎn)進(jìn)行分析前，利用下一代生成矩陣法求出該模型的基本再生數(shù)為

1.2.1 平衡點(diǎn)的非負(fù)性分析

定理1.2.2對(duì)于系統(tǒng)1來說，地方病平衡點(diǎn)E*有以下結(jié)論，當(dāng)R0＞1時(shí)，地方病平衡點(diǎn)E*=（S*，A*，I*，R*）是非負(fù)的。

證明：對(duì)地方病平衡點(diǎn)E*=（S*，A*，I*，R*）進(jìn)行求解，得到

由于參數(shù)非負(fù)，因此S*≥0。對(duì)于A*，I*，R*，根據(jù)等式，可以得知主要取決于λ－dS*。當(dāng)R0≥1時(shí)，同時(shí)因?yàn)?/p>

為證明分?jǐn)?shù)階模型平衡點(diǎn)的穩(wěn)定性可使用如下引理：

引理1[2]分?jǐn)?shù)階系統(tǒng)的平衡點(diǎn)滿足局部漸進(jìn)穩(wěn)定的性質(zhì)，當(dāng)且僅當(dāng)在相應(yīng)平衡點(diǎn)處的Jacobi矩陣J=的所有特征值λi滿足條件。

1.2.2 無病平衡點(diǎn)的局部漸進(jìn)穩(wěn)定性

設(shè)特征值為sα，則模型1（公式（1）所示，下同）在無病平衡點(diǎn)E0處的特征行列式為

定理1.2.3模型1滿足下列結(jié)論，當(dāng)τ≥0時(shí)，該模型的無病平衡點(diǎn)E0有如下結(jié)論：

如果R0＜1，則E0是局部漸進(jìn)穩(wěn)定的。如果R0≥1，則E0是不穩(wěn)定的。

證明：根據(jù)特征矩陣（3），進(jìn)行推導(dǎo)化簡可以得到

根據(jù)（sα+d）2，得到特征值因此。

將行列式部分用下列形式進(jìn)行表示

P1（s）=0存在2個(gè)根，一個(gè)根為一定有負(fù)實(shí)部；另一個(gè)根，當(dāng)R0＜1時(shí)因此此時(shí)它也一定有負(fù)實(shí)部。

1.2.3 地方病平衡點(diǎn)的局部漸進(jìn)穩(wěn)定性

當(dāng)R0＞1時(shí)，模型1存在一個(gè)非負(fù)的地方病平衡點(diǎn)E*=（S*，A*，I*，R*），設(shè)特征值為sα，則地方病平衡點(diǎn)E*處的特征行列式為

先討論τ=0時(shí)，令γ=sα，得到

其中參數(shù)d1，d2，d3可以表示為

而m11，m12，m13，m21，m22，m32，m33可以表示為

因此當(dāng)τ=0時(shí)，E*有如下性質(zhì)：當(dāng)R0＞1時(shí)，如果式（8）滿足n=3時(shí)的霍爾維茨準(zhǔn)則，即式（8）的所有根都具有負(fù)實(shí)部，則E*是局部漸進(jìn)穩(wěn)定的。

當(dāng)τ≠0時(shí)，根據(jù)特征行列式（6），可以化解得到如下的等式

將虛部和實(shí)部分開得到

其中參數(shù)為

對(duì)式（10）經(jīng)過推導(dǎo)化簡，令w2α=y，得到了如下的結(jié)果

進(jìn)一步簡化為

其中各參數(shù)表示意義如下

因此當(dāng)式（13）滿足霍爾維茨準(zhǔn)則時(shí)，式（13）的根均為負(fù)值，則式（10）的解中存在負(fù)實(shí)部，此時(shí)地方病平衡點(diǎn)是局部漸進(jìn)穩(wěn)定的，綜上給出如下定理。

定理1.2.4對(duì)于模型1的地方病平衡點(diǎn)E*，在R0＞1的情況下，若式（8）與式（13）均滿足引理1，則可以得出結(jié)論，當(dāng)τ≥0時(shí)，E*是局部漸進(jìn)穩(wěn)定的。

2 數(shù)值模擬與結(jié)果

為便于后續(xù)計(jì)算，首先通過NSFD算法對(duì)分?jǐn)?shù)階模型進(jìn)行離散化處理[3]，然后進(jìn)行數(shù)據(jù)模擬，主要模擬各類人群隨時(shí)間變化的發(fā)展趨勢。以此來研究帶有時(shí)滯的分?jǐn)?shù)階模型的平衡點(diǎn)的穩(wěn)定性。分2種情況對(duì)模型進(jìn)行數(shù)值模擬，一是對(duì)模型（1）在不同的階數(shù)下隨時(shí)間變化的情況進(jìn)行模擬，二是對(duì)模型在同一分?jǐn)?shù)階下，取不同時(shí)滯時(shí)，隨時(shí)間變化的情況進(jìn)行模擬。

2.1 對(duì)無病平衡點(diǎn)的數(shù)值模擬

由于對(duì)無病平衡點(diǎn)進(jìn)行分析時(shí)要使得不存在非負(fù)的地方病平衡點(diǎn)，因此參數(shù)設(shè)置為λ=5，d=0.04，為增加結(jié)果的可靠性，同時(shí)設(shè)置了幾組不同的階數(shù)進(jìn)行對(duì)比試驗(yàn)，將分?jǐn)?shù)階階數(shù)分別設(shè)置為α=0.7，0.8，0.9，其余參數(shù)分別為β1=β2=0.001，d1=0.03，d2=0.1，δA=0.4，δI=0.3，θ=0.6，τ＝10。根據(jù)參數(shù)設(shè)置，此時(shí)基本再生數(shù)R0≈0.3＜1，且λ-dS*＜0，因此顯然不存在非負(fù)的地方病平衡點(diǎn)，且E0=（125，0，0，0），根據(jù)定理1.2.3，在滿足τ≥10的情況下，當(dāng)R0＜1時(shí)，無病平衡點(diǎn)是局部漸進(jìn)穩(wěn)定的。利用MATLAB進(jìn)行數(shù)值模擬，結(jié)果顯示與理論相符合。因此按照這個(gè)趨勢發(fā)展，各部分人群都將趨于穩(wěn)定，其中只有易感人群還是非零的，其余的都將趨于零。

其余參數(shù)不變，將階數(shù)設(shè)置為0.9，將時(shí)滯分別設(shè)定為τ=5，10，30，此時(shí)模擬結(jié)果各數(shù)值都趨于穩(wěn)定，與定理1.2.3相符合。

2.2 對(duì)地方病平衡點(diǎn)的數(shù)值模擬

在對(duì)地方病平衡點(diǎn)進(jìn)行模擬時(shí)，需要對(duì)參數(shù)重新進(jìn)行修正，選擇參數(shù)λ=10，d=0.02，分?jǐn)?shù)階階數(shù)分別為α=0.7，0.8，0.9，其余參數(shù)保持不變，最終數(shù)值模擬的結(jié)果如圖1所示。

圖1 地方病平衡點(diǎn)下模型各變量在取不同階數(shù)時(shí)隨時(shí)間變化的圖

此時(shí)基本再生數(shù)R0=1.16＞1，顯然存在非負(fù)的地方病平衡點(diǎn)E*=（172，15，23，651）。將假設(shè)參數(shù)帶入式（13）后，根據(jù)定理1.2.4，此時(shí)地方病平衡點(diǎn)是局部漸進(jìn)穩(wěn)定的。顯然數(shù)值模擬的結(jié)果與理論相符合，在該情況下，隨著時(shí)間的發(fā)展，各部分人群都將趨于穩(wěn)定，但是依然存在患者和無癥狀的潛伏者，此時(shí)疫情防控依然不能松懈。

其余參數(shù)不變，將階數(shù)設(shè)置為0.9，將時(shí)滯分別設(shè)定為τ=5，10，30，數(shù)值模擬結(jié)果如圖2所示，此時(shí)，τ≥0，地方病平衡點(diǎn)都是局部漸進(jìn)穩(wěn)定的。模擬的結(jié)果與理論相符合。

圖2 地方病平衡點(diǎn)下模型各變量在取不同時(shí)滯時(shí)隨時(shí)間變化的圖

3 結(jié)論

根據(jù)構(gòu)建的分?jǐn)?shù)階模型，得出了其有2個(gè)平衡點(diǎn)，對(duì)2個(gè)平衡點(diǎn)的局部漸進(jìn)穩(wěn)定性分別進(jìn)行了討論。通過計(jì)算求出了基本再生數(shù)R0，證明了在τ≥0的情況下，無病平衡點(diǎn)E0在R0＜1時(shí)，是局部漸進(jìn)穩(wěn)定的，而當(dāng)R0≥1時(shí)，無病平衡點(diǎn)E0不是局部漸進(jìn)穩(wěn)定的；對(duì)于地方病平衡點(diǎn)E*來說，當(dāng)R0＞1時(shí)，若式（8）與式（13）中參數(shù)均滿足引理1，可以得出結(jié)論，當(dāng)τ≥0時(shí)，E*是局部漸進(jìn)穩(wěn)定的。同時(shí)根據(jù)2種情況對(duì)模型進(jìn)行了模擬，一種是分?jǐn)?shù)階階數(shù)不同，另一種是保持分?jǐn)?shù)階階數(shù)相同但取不同的時(shí)滯，通過模擬驗(yàn)證了之前推導(dǎo)的平衡點(diǎn)穩(wěn)定性的結(jié)論是正確的。此外，通過分析模擬，證明了基本再生數(shù)R0這一指標(biāo)的重要性，為了更好地防控疫情，需要采取措施力求將R0保持在1以下。