基于COVID-19傳染病SIR模型的穩(wěn)定性分析

成紅勝, 成 誠

(1.鹽城師范學(xué)院 信息工程學(xué)院， 江蘇 鹽城 224001; 2.江蘇師范大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院， 江蘇 徐州 221116)

0 引言

2020年年初，一種類似于“非典”的新型冠狀病毒肺炎在武漢被發(fā)現(xiàn)并迅速傳播開來，這場疫情給全人類的生命健康造成了重大威脅.在全民抗擊疫情的過程中，運用數(shù)學(xué)模型預(yù)測疫情的走勢在防控疫情過程中發(fā)揮了重要的作用，為一線的醫(yī)護人員提供了科學(xué)的理論指導(dǎo)和精準的預(yù)測.新冠疫情爆發(fā)已經(jīng)過去一年多，世界各地陸續(xù)發(fā)現(xiàn)感染新冠病毒康復(fù)后復(fù)陽的病例.基于這種實際情況，研究一類具有飽和發(fā)生率和暫時免疫力的傳染病SIR模型，得到平衡點的相關(guān)性質(zhì)，有望從科學(xué)的角度研究傳染病的傳播機制.

長久以來，人們對于傳染病模型的研究已經(jīng)取得了一些顯著的成果.Wen等學(xué)者研究了一類具有暫時免疫力的SIR模型[1]：

(1)

其中:S(t)的表示容易被傳染該種疾病的總?cè)藬?shù)；I(t)表示已經(jīng)被傳染該種疾病的人數(shù)；R(t)表示因接種疫苗從而獲得暫時免疫力的人數(shù);μ表示由易感染變成感染者的速率.參數(shù)μ1、μ2、μ3均為非負常數(shù)，μ1表示易感染人群的死亡率；μ2表示感染者的死亡率；μ3表示從該種傳染病康復(fù)者的死亡率.從生物學(xué)的角度來看，μ1

(2)

系統(tǒng)(2)中的初始條件表示如下：

(3)

1 無病平衡點的穩(wěn)定性

(4)

分離變量可得

(5)

等式兩邊同時積分可得方程的解為

(6)

假設(shè)ξ表示感染者的平均壽命，可得概率函數(shù)為

(7)

其中的概率密度函數(shù)為

(8)

由此計算得到數(shù)學(xué)期望

(9)

因此系統(tǒng)(2)存在一個無病平衡點E0.通過計算E0處的Jacobi矩陣，得到特征方程為

(10)

求解得到三個特征值分別為

定理3 當(dāng)R0<1且μ1=μ2=μ3時，系統(tǒng)(2)在E0(S0,0,0)處在全局的性質(zhì)是漸近穩(wěn)定的.

證明構(gòu)造如下的Lyapunov函數(shù)：

(11)

等式兩端同時對時間t求導(dǎo)得

[S(t)+I(t)+R(t)-S0](μ-μ1S(t)-μ2I(t)-μ3R(t))

(12)

因為μ1=μ2=μ3，進一步簡化式(12)得

[S(t)+I(t)+R(t)-S0](μS(0)-μ1S(t)-μ2I(t)=

-μ1[S(t)+I(t)+R(t)-S0]2=

-μ1[[S(t)-S0]+I(t)+R(t)]2

(13)

2 地方病平衡點的穩(wěn)定性

定理4 系統(tǒng)(2)存在一個地方病平衡點E*(S*,I*,R*)，當(dāng)R0>1時，系統(tǒng)(2)在E*處是局部性質(zhì)漸近穩(wěn)定的；當(dāng)R0<1時，系統(tǒng)(2)在E*處是不穩(wěn)定的.

證明證明方法與定理2的證明過程類似.

令系統(tǒng)(2)中等式左側(cè)等于零，即可求得這個系統(tǒng)的地方病平衡點為E*(S*,I*,R*)， 具體為

通過計算系統(tǒng)(2)在E*處的Jacobi矩陣，得到特征方程滿足如下的形式：

[λ2+p0(τ)+p1(τ)λ+q0(τ)eλτ](λ+μ3)=0

(14)

其中

顯然，特征方程(14)有一個特征值-μ3且其實部為負數(shù). 其余的特征值由以下方程所確定

λ2+p0(τ)+p1(τ)λ+q0(τ)eλτ

(15)

下面將τ分為τ=0和τ>0兩種情況進行討論：

1) 當(dāng)τ=0時，方程(15)可以寫成：

λ2+p1(0)λ+p0(0)+q0(0)=0

(16)

2) 當(dāng)τ>0時，假設(shè)iω(ω>0)是方程(15)的一個解，將iω代入方程，分離實部和虛部得到如下的等價方程：

(17)

將兩個方程分別平方再相加，得到

(18)

其中

因此

綜合1)和2)兩種情況可知，當(dāng)R0>1，對任意的τ≥0，地方病平衡點E*(S*,I*,R*)是局部漸近穩(wěn)定的.當(dāng)R0<1時，特征方程(15)必有一解是正數(shù)，故存在一個特征值的實部是正數(shù)，此時系統(tǒng)在E*處是不穩(wěn)定的.證畢.

定理5當(dāng)R0>1且μ1=μ2=μ3時，系統(tǒng)(2)在E*(S*,I*.R*)處是全局漸近穩(wěn)定的.

證明構(gòu)造Lyapunov泛函：

(19)

等式兩端同時對時間t求導(dǎo)得

[(S(t)-S*)+(I(t)-I*)+(R(t)-R*)](μ-μ1S(t)-βS(t)I(t)+

γI(t-τ)e-μ3τ+βS(t)I(t)-(μ2+γ)I(t)+γI(t)-μ3R(t)-γI(t-τ)e-μ3τ=

[(S(t)-S*)+(I(t)-I*)+(R(t)-R*)](μ-μ1S(t)-μ2I(t)-μ3R(t))

將μ1=μ2=μ3代入，上面的式子可以寫成

[(S(t)-S*)+(I(t)-I*)+(R(t)-R*)][μ-μ(S(t)-S*)-μ(I(t)-I*)-

μ(R(t)-R*)-μ(S*+I*+R*)]=

[(S(t)-S*)+(I(t)-I*)+(R(t)-R*)][-μ(S(t)-S*)+(I(t)-I*)+

(R(t)-R*)]=-μ[(S(t)-S*)+(I(t)-I*)+(R(t)-R*)]2

3 全局穩(wěn)定性的數(shù)值模擬

1) 選取μ=0.3，α=2，β=0.1，γ=0.3，μ1=0.1，μ2=0.1，μ3=0.1，τ=15.計算得到R0=0.6<1，無病平衡點E0=(3,0,0).由定理3可知，系統(tǒng)(2)在E0處是全局漸近穩(wěn)定的.通過Matlab進行數(shù)值模擬驗證了該結(jié)論，結(jié)果如圖1所示.

圖1 無病平衡點的穩(wěn)定性 圖2 地方病平衡點的穩(wěn)定性

4 結(jié)論

基于疫情常態(tài)化的實際情況，研究了一種傳染病SIR模型，這種模型具有飽和發(fā)生率和暫時免疫力. 通過分析特征函數(shù)和特征值，得到當(dāng)R0<1時，無病平衡點E0在局部的性質(zhì)是漸近穩(wěn)定的，該傳染病并不會傳播，也不會變成一種流行??；當(dāng)R0>1時，地方病平衡點E*在局部的性質(zhì)是漸近穩(wěn)定的，這種傳染病最終會傳播開來.通過建立Lyapunov函數(shù)，得到當(dāng)R0<1且μ1=μ2=μ3時，無病平衡點E0在全局的性質(zhì)是漸近穩(wěn)定的；當(dāng)R0>1且μ1=μ2=μ3時，地方病平衡點在全局的性質(zhì)是漸近穩(wěn)定的.最后通過Matlab進行模擬，驗證所得到的結(jié)果.