COVID-19分數(shù)階群體時滯模型的穩(wěn)定性分析

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

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

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

1 模型的建立與分析

1.1 分數(shù)階模型的建立

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

上述系統(tǒng)當t=0時，會滿足下述的初值條件

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

1.2 模型的分析

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

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

1.2.1 平衡點的非負性分析

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

證明：對地方病平衡點E*=（S*，A*，I*，R*）進行求解，得到

由于參數(shù)非負，因此S*≥0。對于A*，I*，R*，根據(jù)等式，可以得知主要取決于λ－dS*。當R0≥1時，同時因為

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

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

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

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

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

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

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

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

將行列式部分用下列形式進行表示

P1（s）=0存在2個根，一個根為一定有負實部；另一個根，當R0＜1時因此此時它也一定有負實部。

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

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

先討論τ=0時，令γ=sα，得到

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

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

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

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

將虛部和實部分開得到

其中參數(shù)為

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

進一步簡化為

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

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

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

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

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

2.1 對無病平衡點的數(shù)值模擬

由于對無病平衡點進行分析時要使得不存在非負的地方病平衡點，因此參數(shù)設(shè)置為λ=5，d=0.04，為增加結(jié)果的可靠性，同時設(shè)置了幾組不同的階數(shù)進行對比試驗，將分數(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ù)R0≈0.3＜1，且λ-dS*＜0，因此顯然不存在非負的地方病平衡點，且E0=（125，0，0，0），根據(jù)定理1.2.3，在滿足τ≥10的情況下，當R0＜1時，無病平衡點是局部漸進穩(wěn)定的。利用MATLAB進行數(shù)值模擬，結(jié)果顯示與理論相符合。因此按照這個趨勢發(fā)展，各部分人群都將趨于穩(wěn)定，其中只有易感人群還是非零的，其余的都將趨于零。

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

2.2 對地方病平衡點的數(shù)值模擬

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

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

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

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

圖2 地方病平衡點下模型各變量在取不同時滯時隨時間變化的圖

3 結(jié)論

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